###
# File description: RICLPM and CLPM analysis
# Script author: Justin Robinson
# Last modified: 12-01-2024
###

#### 1. Set-up and data loading ----------

rm( list = ls() ) # clearing our environment 

## Loading packages 

library(lavaan)
library(OrthoPanels)
library(semTools)
library(foreign)
library(dplyr)
library(readstata13)
library(rstudioapi)   # for automatic setwd()

setwd(dirname(getActiveDocumentContext()$path))

getwd()


dat <- read.csv("RICLPM_data.csv")


#### 2. Measurement models ---------

# Authoritarianism

auth_MM <- 
  "
auth7 =~ auth1_W7 + auth2_W7 + auth3_W7 + auth4_W7
auth10 =~ auth1_W10 + auth2_W10 + auth3_W10 + auth4_W10
auth11 =~ auth1_W11 + auth2_W11 + auth3_W11 + auth4_W11
auth14 =~ auth1_W14 + auth2_W14 + auth3_W14 + auth4_W14
auth10 ~ auth7
auth11 ~ auth10
auth14 ~ auth11
"
auth_MM_fit <- sem(auth_MM, 
                   test="Satorra-Bentler",
                   ordered = c("auth1_W7",
                               "auth2_W7",
                               "auth3_W7",
                               "auth4_W7",
                               "auth1_W10",
                               "auth2_W10",
                               "auth3_W10",
                               "auth4_W10",
                               "auth1_W11",
                               "auth2_W11",
                               "auth3_W11",
                               "auth4_W11",
                               "auth1_W14",
                               "auth2_W14",
                               "auth3_W14",
                               "auth4_W14"),
                   data=dat)
reliability(auth_MM_fit, return.total = TRUE)
fitMeasures(auth_MM_fit, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
auth_MM_results <- parameterEstimates(auth_MM_fit)


auth_MM_2 <- 
  "
auth7 =~ auth1_W7 + auth2_W7 + auth3_W7 + auth4_W7
auth10 =~ auth1_W10 + auth2_W10 + auth3_W10 + auth4_W10
auth11 =~ auth1_W11 + auth2_W11 + auth3_W11 + auth4_W11
auth14 =~ auth1_W14 + auth2_W14 + auth3_W14 + auth4_W14
auth1_W7 ~~ auth1_W10 + auth1_W11 + auth1_W14
auth1_W10 ~~ auth1_W11 + auth1_W14
auth1_W11 ~~ auth1_W14
auth2_W7 ~~ auth2_W10 + auth2_W11 + auth2_W14
auth2_W10 ~~ auth2_W11 + auth2_W14
auth2_W11 ~~ auth2_W14
auth3_W7 ~~ auth3_W10 + auth3_W11 + auth3_W14
auth3_W10 ~~ auth3_W11 + auth3_W14
auth3_W11 ~~ auth3_W14
auth4_W7 ~~ auth4_W10 + auth4_W11 + auth4_W14
auth4_W10 ~~ auth4_W11 + auth4_W14
auth4_W11 ~~ auth4_W14

auth10 ~ auth7
auth11 ~ auth10
auth14 ~ auth11

"
auth_MM_2_fit <- sem(auth_MM_2, 
                     test="Satorra-Bentler",
                     ordered = c("auth1_W7",
                                 "auth2_W7",
                                 "auth3_W7",
                                 "auth4_W7",
                                 "auth1_W10",
                                 "auth2_W10",
                                 "auth3_W10",
                                 "auth4_W10",
                                 "auth1_W11",
                                 "auth2_W11",
                                 "auth3_W11",
                                 "auth4_W11",
                                 "auth1_W14",
                                 "auth2_W14",
                                 "auth3_W14",
                                 "auth4_W14"),
                     data=dat)
reliability(auth_MM_2_fit, return.total = TRUE)
summary(auth_MM_2_fit)
fitMeasures(auth_MM_2_fit, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))
auth_MM_2_results <- parameterEstimates(auth_MM_2_fit)

anova(auth_MM_fit, auth_MM_2_fit)

## Immigratio

imm_pref_2_MM <- 
  "
imm7 =~ im_culturalW7 + immigSelf01W7
imm10 =~ im_culturalW10 + immigSelf01W10
imm11 =~ im_culturalW11 + immigSelf01W11
imm14 =~ im_culturalW14 + immigSelf01W14
imm10 ~ imm7
imm11 ~ imm10
imm14 ~ imm11
"
imm_pref_2_MM_fit <- sem(imm_pref_2_MM, 
                         test="Satorra-Bentler",
                         missing = "ML",
                         data=dat)
reliability(imm_pref_2_MM_fit, return.total = TRUE)
summary(imm_pref_2_MM_fit)
fitMeasures(imm_pref_2_MM_fit, c("chisq.scaled", "df", "pvalue.scaled", "cfi", "tli", "rmsea", "srmr", "aic"))

## Redistribution

redist_MM <- 
  "
red7 =~ 1*redistSelf01W7
red10 =~ 1*redistSelf01W10
red11 =~ 1*redistSelf01W11
red14 =~ 1*redistSelf01W14
red10 ~ red7
red11 ~ red10
red14 ~ red11
"
redist_MM_fit <- sem(redist_MM, se="robust.huber.white", test="Satorra-Bentler", data=dat)
reliability(redist_MM_fit, treturn.total = TRUE)
reliabilityL2(redist_MM_fit, "red7")
reliabilityL2(redist_MM_fit, "red10")
reliabilityL2(redist_MM_fit, "red11")
reliabilityL2(redist_MM_fit, "red14")

summary(redist_MM_fit)

## EU integration

eu_int_MM <- "
EU7 =~ eu_intW7
EU10 =~ eu_intW10
EU11 =~ eu_intW11
EU14 =~ eu_intW14
EU10 ~ EU7
EU11 ~ EU10
EU14 ~ EU11
"
eu_int_MM_fit <- sem(eu_int_MM, 
                     test="Satorra-Bentler", 
                     data=dat)
reliability(eu_int_MM_fit, return.total = TRUE)
reliabilityL2(eu_int_MM_fit, "EU7")
reliabilityL2(eu_int_MM_fit, "EU10")
reliabilityL2(eu_int_MM_fit, "EU11")
reliabilityL2(eu_int_MM_fit, "EU14")
summary(eu_int_MM_fit)

#### 3. Measurement invariance ---------

## Fitting the models 

configural <- '
  RIauth =~ 1*auth01W7 + 1*auth01W10 + 1*auth01W11 + 1*auth01W14
  RIim_prefW =~ 1*im_pref_2W7 + 1*im_pref_2W10 + 1*im_pref_2W11 + 1*im_pref_2W14  

  # Create within-person parts
  wauth01W7 =~ a*auth01W7
  wauth01W10 =~ a*auth01W10
  wauth01W11 =~ a*auth01W11 
  wauth01W14 =~ a*auth01W14
  
  wim_prefW7 =~ b*im_pref_2W7
  wim_prefW10 =~ b*im_pref_2W10
  wim_prefW11 =~ b*im_pref_2W11
  wim_prefW14 =~ b*im_pref_2W14
  
   # Estimate the lagged effects between the within-person centered variables.
  wauth01W10 + wim_prefW10 ~ wauth01W7 + wim_prefW7
  wauth01W11 + wim_prefW11 ~ wauth01W10 + wim_prefW10
  wauth01W14 + wim_prefW14 ~ wauth01W11 + wim_prefW11

  # Estimate the covariance between the within-person centered variables
  wauth01W7 ~~ wim_prefW7 # Covariance
  wauth01W10 ~~ wim_prefW10
  wauth01W11 ~~ wim_prefW11
  wauth01W14 ~~ wim_prefW14
  
# Constrain cov of between factors and exogenous within factors to 0:
  RIauth + RIim_prefW ~~ 0*wauth01W7 + 0*wim_prefW7
  
    # Fix error variance of indicators to (1 - alpha)*var
auth01W7 ~~ 0.02*auth01W7
auth01W10 ~~ 0.016*auth01W10
auth01W11 ~~ 0.018*auth01W11
auth01W14 ~~ 0.019*auth01W14
im_pref_2W7 ~~ 0.0147*im_pref_2W7
im_pref_2W10 ~~ 0.000155*im_pref_2W10
im_pref_2W11 ~~ 0.000147*im_pref_2W11
im_pref_2W14 ~~ 0.000512*im_pref_2W14

  '
configural.fit <- cfa(configural, 
                      data = dat, 
                      missing = 'ML',
                      estimator = "MLR")
summary(configural.fit, standardized = T)


weak_fa <- '
  RIauth =~ 1*auth01W7 + 1*auth01W10 + 1*auth01W11 + 1*auth01W14
  RIim_prefW =~ 1*im_pref_2W7 + 1*im_pref_2W10 + 1*im_pref_2W11 + 1*im_pref_2W14  

  # Create within-person parts
  wauth01W7 =~ 1*auth01W7
  wauth01W10 =~ 1*auth01W10
  wauth01W11 =~ 1*auth01W11 
  wauth01W14 =~ 1*auth01W14
  
  wim_prefW7 =~ 1*im_pref_2W7
  wim_prefW10 =~ 1*im_pref_2W10
  wim_prefW11 =~ 1*im_pref_2W11
  wim_prefW14 =~ 1*im_pref_2W14
  
   # Estimate the lagged effects between the within-person centered variables.
  wauth01W10 + wim_prefW10 ~ wauth01W7 + wim_prefW7
  wauth01W11 + wim_prefW11 ~ wauth01W10 + wim_prefW10
  wauth01W14 + wim_prefW14 ~ wauth01W11 + wim_prefW11

  # Estimate the covariance between the within-person centered variables
  wauth01W7 ~~ wim_prefW7 # Covariance
  wauth01W10 ~~ wim_prefW10
  wauth01W11 ~~ wim_prefW11
  wauth01W14 ~~ wim_prefW14
  
# Constrain cov of between factors and exogenous within factors to 0:
  RIauth + RIim_prefW ~~ 0*wauth01W7 + 0*wim_prefW7
  
    # Fix error variance of indicators to (1 - alpha)*var
auth01W7 ~~ 0.02*auth01W7
auth01W10 ~~ 0.016*auth01W10
auth01W11 ~~ 0.018*auth01W11
auth01W14 ~~ 0.019*auth01W14
im_pref_2W7 ~~ 0.0147*im_pref_2W7
im_pref_2W10 ~~ 0.000155*im_pref_2W10
im_pref_2W11 ~~ 0.000147*im_pref_2W11
im_pref_2W14 ~~ 0.000512*im_pref_2W14

  '
weak_fa.fit <- cfa(weak_fa, 
                   data = dat, 
                   missing = 'ML',
                   estimator = "MLR")
summary(weak_fa.fit, standardized = T)

strong_fa <- '
  RIauth =~ 1*auth01W7 + 1*auth01W10 + 1*auth01W11 + 1*auth01W14
  RIim_prefW =~ 1*im_pref_2W7 + 1*im_pref_2W10 + 1*im_pref_2W11 + 1*im_pref_2W14  

  # Create within-person parts
  wauth01W7 =~ 1*auth01W7
  wauth01W10 =~ 1*auth01W10
  wauth01W11 =~ 1*auth01W11 
  wauth01W14 =~ 1*auth01W14
  
  wim_prefW7 =~ 1*im_pref_2W7
  wim_prefW10 =~ 1*im_pref_2W10
  wim_prefW11 =~ 1*im_pref_2W11
  wim_prefW14 =~ 1*im_pref_2W14
  
  # Constrained intercepts over time
  auth01W7 + auth01W10 + auth01W11 + auth01W14 ~ c*1
  im_pref_2W7 + im_pref_2W10 + im_pref_2W11 + im_pref_2W14 ~ d*1 
  
  # Free latent means from t = 2 onward
  wauth01W10 + wauth01W11 + wauth01W14 + wim_prefW10 + wim_prefW11 + wim_prefW14 ~ 1
  
   # Estimate the lagged effects between the within-person centered variables.
  wauth01W10 + wim_prefW10 ~ wauth01W7 + wim_prefW7
  wauth01W11 + wim_prefW11 ~ wauth01W10 + wim_prefW10
  wauth01W14 + wim_prefW14 ~ wauth01W11 + wim_prefW11

  # Estimate the covariance between the within-person centered variables
  wauth01W7 ~~ wim_prefW7 # Covariance
  wauth01W10 ~~ wim_prefW10
  wauth01W11 ~~ wim_prefW11
  wauth01W14 ~~ wim_prefW14
  
# Constrain cov of between factors and exogenous within factors to 0:
  RIauth + RIim_prefW ~~ 0*wauth01W7 + 0*wim_prefW7
  
    # Fix error variance of indicators to (1 - alpha)*var
auth01W7 ~~ 0.02*auth01W7
auth01W10 ~~ 0.016*auth01W10
auth01W11 ~~ 0.018*auth01W11
auth01W14 ~~ 0.019*auth01W14
im_pref_2W7 ~~ 0.0147*im_pref_2W7
im_pref_2W10 ~~ 0.000155*im_pref_2W10
im_pref_2W11 ~~ 0.000147*im_pref_2W11
im_pref_2W14 ~~ 0.000512*im_pref_2W14

  '
strong_fa.fit <- cfa(strong_fa, 
                     data = dat, 
                     missing = 'ML',
                     estimator = "MLR")
summary(strong_fa.fit, standardized = T)

# Comparing model fit
fitmeasures(configural.fit)
fitmeasures(weak_fa.fit)
fitmeasures(strong_fa.fit)

#### 4. CL2PMs ----------------

## Immigration

immig_pref_2_CL2PM <- '
   # Create between components (random intercepts)
  RIauth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wauth01W7 =~ 1*std_auth01W7
  wauth01W10 =~ 1*std_auth01W10
  wauth01W11 =~ 1*std_auth01W11 
  wauth01W14 =~ 1*std_auth01W14

  wim_prefW7 =~ 1*std_im_pref_2W7
  wim_prefW10 =~ 1*std_im_pref_2W10
  wim_prefW11 =~ 1*std_im_pref_2W11
  wim_prefW14 =~ 1*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wauth01W10 + wim_prefW10 ~ wauth01W7 + wim_prefW7
  wauth01W11 + wim_prefW11 ~ wauth01W10 + wim_prefW10 + wauth01W7 + wim_prefW7
  wauth01W14 + wim_prefW14 ~ wauth01W11 + wim_prefW11 + wauth01W10 + wim_prefW10

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wauth01W7 ~~ wim_prefW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wauth01W10 ~~ wim_prefW10
  wauth01W11 ~~ wim_prefW11
  wauth01W14 ~~ wim_prefW14

  # Estimate the (residual) variance of the within-person centered variables.
  wauth01W7 ~~ wauth01W7 # Variances
  wim_prefW7 ~~ wim_prefW7 
  wauth01W10 ~~ wauth01W10 # Residual variances
  wim_prefW10 ~~ wim_prefW10 
  wauth01W11 ~~ wauth01W11
  wim_prefW11 ~~ wim_prefW11 
  wauth01W14 ~~ wauth01W14 
  wim_prefW14 ~~ wim_prefW14 
  
  # Fix the variance and covariance of the random intercepts to 0:
  RIauth01W ~~ 0*RIauth01W
  RIim_prefW ~~ 0*RIim_prefW
  RIauth01W ~~ 0*RIim_prefW
  
  '
immig_pref_2_CL2PM.fit <- lavaan(immig_pref_2_CL2PM, data = dat, 
                                 estimator = "MLR",
                                 missing = 'ML', meanstructure = T, int.ov.free = T) 
summary(immig_pref_2_CL2PM.fit, standardized = T)

## EU integration

std_eu_pref_CL2PM <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wstd_eu_intW7 =~ 1*std_eu_intW7
  wstd_eu_intW10 =~ 1*std_eu_intW10
  wstd_eu_intW11 =~ 1*std_eu_intW11
  wstd_eu_intW14 =~ 1*std_eu_intW14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 + wstd_eu_intW10 ~ wstd_auth01W7 + wstd_eu_intW7
  wstd_auth01W11 + wstd_eu_intW11 ~ wstd_auth01W10 + wstd_eu_intW10 + wstd_auth01W7 + wstd_eu_intW7
  wstd_auth01W14 + wstd_eu_intW14 ~ wstd_auth01W11 + wstd_eu_intW11 + wstd_auth01W10 + wstd_eu_intW10

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wstd_eu_intW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wstd_eu_intW10
  wstd_auth01W11 ~~ wstd_eu_intW11
  wstd_auth01W14 ~~ wstd_eu_intW14
  
  
  # Constrain the variance and covariance of the random intercepts to zero 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIstd_eu_intW ~~ 0*RIstd_eu_intW
  RIstd_auth01W ~~ 0*RIstd_eu_intW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ wstd_eu_intW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ wstd_eu_intW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wstd_eu_intW11 ~~ wstd_eu_intW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wstd_eu_intW14 ~~ wstd_eu_intW14 
  

'

std_eu_pref_CL2PM.fit <- lavaan(std_eu_pref_CL2PM, data = dat, 
                                estimator = "MLR",
                                missing = 'ML', meanstructure = T, int.ov.free = T) 
summary(std_eu_pref_CL2PM.fit, standardized = T)

standardizedsolution(std_eu_pref_CL2PM.fit)

#### 5. RI-CLPMs ----------

# I fit the models here, but do not report the results from these models. 
# Instead, model results are derived from the standardized solution models 
# which are presented in section X. 


### Immigration

std_immig_pref_2_tvc <- '
  # Create between components (random intercepts)
  RIauth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wauth01W7 =~ 1*std_auth01W7
  wauth01W10 =~ 1*std_auth01W10
  wauth01W11 =~ 1*std_auth01W11 
  wauth01W14 =~ 1*std_auth01W14

  wim_prefW7 =~ 1*std_im_pref_2W7
  wim_prefW10 =~ 1*std_im_pref_2W10
  wim_prefW11 =~ 1*std_im_pref_2W11
  wim_prefW14 =~ 1*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wauth01W10 + wim_prefW10 ~ wauth01W7 + wim_prefW7
  wauth01W11 + wim_prefW11 ~ wauth01W10 + wim_prefW10
  wauth01W14 + wim_prefW14 ~ wauth01W11 + wim_prefW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wauth01W7 ~~ wim_prefW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wauth01W10 ~~ wim_prefW10
  wauth01W11 ~~ wim_prefW11
  wauth01W14 ~~ wim_prefW14
  
  
  # Estimate the variance and covariance of the random intercepts. 
  RIauth01W ~~ RIauth01W
  RIim_prefW ~~ RIim_prefW
  RIauth01W ~~ RIim_prefW

  # Estimate the (residual) variance of the within-person centered variables.
  wauth01W7 ~~ wauth01W7 # Variances
  wim_prefW7 ~~ wim_prefW7 
  wauth01W10 ~~ wauth01W10 # Residual variances
  wim_prefW10 ~~ wim_prefW10 
  wauth01W11 ~~ wauth01W11
  wim_prefW11 ~~ wim_prefW11 
  wauth01W14 ~~ wauth01W14 
  wim_prefW14 ~~ wim_prefW14 
  
# Fix error variance of indicators to (1 - alpha)*var
std_auth01W7 ~~ 0.02*std_auth01W7
std_auth01W10 ~~ 0.016*std_auth01W10
std_auth01W11 ~~ 0.018*std_auth01W11
std_auth01W14 ~~ 0.019*std_auth01W14
std_im_pref_2W7 ~~ 0.0147*std_im_pref_2W7
std_im_pref_2W10 ~~ 0.000155*std_im_pref_2W10
std_im_pref_2W11 ~~ 0.000147*std_im_pref_2W11
std_im_pref_2W14 ~~ 0.000512*std_im_pref_2W14

# incorporate tvc - age
wauth01W7 + wim_prefW7 ~ ageW7
wauth01W10 + wim_prefW10 ~ ageW10
wauth01W11 + wim_prefW11 ~ ageW11
wauth01W14 + wim_prefW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc – social grade
wauth01W7 + wim_prefW7 ~ p_socgradeW7
wauth01W10 + wim_prefW10 ~ p_socgradeW10
wauth01W11 + wim_prefW11 ~ p_socgradeW11
wauth01W14 + wim_prefW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – income
wauth01W7 + wim_prefW7 ~ incomerec7
wauth01W10 + wim_prefW10 ~ incomerec10
wauth01W11 + wim_prefW11 ~ incomerec11
wauth01W14 + wim_prefW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14
  '
std_immig_pref_2_tvc.fit <- lavaan(std_immig_pref_2_tvc, 
                                   data = dat, 
                                   estimator = "MLR",
                                   missing = 'ML', 
                                   meanstructure = T, 
                                   int.ov.free = T) 
summary(std_immig_pref_2_tvc.fit, standardized = T)

# EU integration

std_eu_pref_tvc <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wstd_eu_intW7 =~ 1*std_eu_intW7
  wstd_eu_intW10 =~ 1*std_eu_intW10
  wstd_eu_intW11 =~ 1*std_eu_intW11
  wstd_eu_intW14 =~ 1*std_eu_intW14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 + wstd_eu_intW10 ~ wstd_auth01W7 + wstd_eu_intW7
  wstd_auth01W11 + wstd_eu_intW11 ~ wstd_auth01W10 + wstd_eu_intW10
  wstd_auth01W14 + wstd_eu_intW14 ~ wstd_auth01W11 + wstd_eu_intW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wstd_eu_intW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wstd_eu_intW10
  wstd_auth01W11 ~~ wstd_eu_intW11
  wstd_auth01W14 ~~ wstd_eu_intW14
  
  
  # Estimate the variance and covariance of the random intercepts. 
  RIstd_auth01W ~~ RIstd_auth01W
  RIstd_eu_intW ~~ RIstd_eu_intW
  RIstd_auth01W ~~ RIstd_eu_intW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ wstd_eu_intW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ wstd_eu_intW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wstd_eu_intW11 ~~ wstd_eu_intW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wstd_eu_intW14 ~~ wstd_eu_intW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W10 ~~ 0.0058*std_auth01W10
std_auth01W11 ~~ 0.002*std_auth01W11
std_auth01W14 ~~ 0.00725*std_auth01W14
std_eu_intW10 ~~ 0.0028*std_eu_intW10
std_eu_intW11 ~~ 0.00194*std_eu_intW11
std_eu_intW14 ~~ 0*std_eu_intW14

# incorporate tvc - age
wstd_auth01W7 + wstd_eu_intW7 ~ ageW7
wstd_auth01W10 + wstd_eu_intW10 ~ ageW10
wstd_auth01W11 + wstd_eu_intW11 ~ ageW11
wstd_auth01W14 + wstd_eu_intW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc – income
wstd_auth01W7 + wstd_eu_intW7 ~ incomerec7
wstd_auth01W10 + wstd_eu_intW10 ~ incomerec10
wstd_auth01W11 + wstd_eu_intW11 ~ incomerec11
wstd_auth01W14 + wstd_eu_intW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14

# incorporate tvc – social grade
wstd_auth01W7 + wstd_eu_intW7 ~ p_socgradeW7
wstd_auth01W10 + wstd_eu_intW10 ~ p_socgradeW10
wstd_auth01W11 + wstd_eu_intW11 ~ p_socgradeW11
wstd_auth01W14 + wstd_eu_intW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

'

std_eu_pref_tvc.fit <- lavaan(std_eu_pref_tvc, data = dat, 
                              estimator = "MLR",
                              missing = 'ML', meanstructure = T, int.ov.free = T) 
summary(std_eu_pref_tvc.fit, standardized = T)

#### 6. RI-CLPMs - equality constraints ----------

std_immig_pref_EQ <- '
 # Create between components (random intercepts)
  RIauth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wauth01W7 =~ NA*std_auth01W7
  wauth01W10 =~ NA*std_auth01W10
  wauth01W11 =~ NA*std_auth01W11 
  wauth01W14 =~ NA*std_auth01W14

  wim_prefW7 =~ NA*std_im_pref_2W7
  wim_prefW10 =~ NA*std_im_pref_2W10
  wim_prefW11 =~ NA*std_im_pref_2W11
  wim_prefW14 =~ NA*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wauth01W10 ~ a*wauth01W7 + b*wim_prefW7
  wim_prefW10 ~ c*wauth01W7 + d*wim_prefW7
  wauth01W11 ~ a*wauth01W10 + b*wim_prefW10
  wim_prefW11 ~ c*wauth01W10 + d*wim_prefW10
  wauth01W14 ~ a*wauth01W11 + b*wim_prefW11
  wim_prefW14 ~ c*wauth01W11 + d*wim_prefW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wauth01W7 ~~ cor1*wim_prefW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wauth01W10 ~~ rcov2*wim_prefW10
  wauth01W11 ~~ rcov3*wim_prefW11
  wauth01W14 ~~ rcov4*wim_prefW14
 
  # Estimate the variance and covariance of the random intercepts. 
  RIauth01W ~~ 0.691*RIauth01W
  RIim_prefW ~~ 0.746*RIim_prefW
  RIauth01W ~~ RIim_prefW

  # Estimate the (residual) variance of the within-person centered variables.
  wauth01W7 ~~ 1*wauth01W7 # Variances
  wim_prefW7 ~~ 1*wim_prefW7 
  
  # Label residual variances
  wauth01W10 ~~ rvx2*wauth01W10 # Residual variances
  wim_prefW10 ~~ rvy2*wim_prefW10 
  wauth01W11 ~~ rvx3*wauth01W11
  wim_prefW11 ~~ rvy3*wim_prefW11 
  wauth01W14 ~~ rvx4*wauth01W14 
  wim_prefW14 ~~ rvy4*wim_prefW14 
  
# Fix error variance of indicators to (1 - alpha)*var
std_auth01W7 ~~ 0.02*std_auth01W7
std_auth01W10 ~~ 0.016*std_auth01W10
std_auth01W11 ~~ 0.018*std_auth01W11
std_auth01W14 ~~ 0.019*std_auth01W14
std_im_pref_2W7 ~~ 0.0147*std_im_pref_2W7
std_im_pref_2W10 ~~ 0.000155*std_im_pref_2W10
std_im_pref_2W11 ~~ 0.000147*std_im_pref_2W11
std_im_pref_2W14 ~~ 0.000512*std_im_pref_2W14

# incorporate tvc - age
wauth01W7 + wim_prefW7 ~ ageW7
wauth01W10 + wim_prefW10 ~ ageW10
wauth01W11 + wim_prefW11 ~ ageW11
wauth01W14 + wim_prefW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc – social grade
wauth01W7 + wim_prefW7 ~ p_socgradeW7
wauth01W10 + wim_prefW10 ~ p_socgradeW10
wauth01W11 + wim_prefW11 ~ p_socgradeW11
wauth01W14 + wim_prefW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – income
wauth01W7 + wim_prefW7 ~ incomerec7
wauth01W10 + wim_prefW10 ~ incomerec10
wauth01W11 + wim_prefW11 ~ incomerec11
wauth01W14 + wim_prefW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14

 # Compute correlations of within-components at each wave
  cor2 := a*c + b*d + a*d*cor1 + b*c*cor1 + rcov2
  cor3 := a*c + b*d + a*d*cor2 + b*c*cor2 + rcov3
  cor4 := a*c + b*d + a*d*cor3 + b*c*cor3 + rcov4
 
  # Contrain residual variances of within-components such that variance of each 
  # within-component equals 1
  rvx2 == 1 - (a*a + b*b + 2*a*b*cor1)
  rvy2 == 1 - (c*c + d*d + 2*c*d*cor1)
  rvx3 == 1 - (a*a + b*b + 2*a*b*cor2)
  rvy3 == 1 - (c*c + d*d + 2*c*d*cor2)
  rvx4 == 1 - (a*a + b*b + 2*a*b*cor3)
  rvy4 == 1 - (c*c + d*d + 2*c*d*cor3)
  '

std_immig_pref_EQ.fit <- lavaan(std_immig_pref_EQ, 
                                data = dat, 
                                estimator = "MLR",
                                missing = 'ML', 
                                meanstructure = T, 
                                int.ov.free = T) 
summary(std_immig_pref_EQ.fit, standardized = T)

std_eu_EQ <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ NA*std_auth01W7
  wstd_auth01W10 =~ NA*std_auth01W10
  wstd_auth01W11 =~ NA*std_auth01W11 
  wstd_auth01W14 =~ NA*std_auth01W14

  wstd_eu_intW7 =~ NA*std_eu_intW7
  wstd_eu_intW10 =~ NA*std_eu_intW10
  wstd_eu_intW11 =~ NA*std_eu_intW11
  wstd_eu_intW14 =~ NA*std_eu_intW14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 ~ a*wstd_auth01W7 + b*wstd_eu_intW7
  wstd_eu_intW10 ~ c*wstd_auth01W7 + d*wstd_eu_intW7
  wstd_auth01W11 ~ a*wstd_auth01W10 + b*wstd_eu_intW10
  wstd_eu_intW11 ~ c*wstd_auth01W10 + d*wstd_eu_intW10
  wstd_auth01W14 ~ a*wstd_auth01W11 + b*wstd_eu_intW11
  wstd_eu_intW14 ~ c*wstd_auth01W11 + d*wstd_eu_intW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ cor1*wstd_eu_intW7 # Covariance
  
  # Label residual covariances
  wstd_auth01W10 ~~ rcov2*wstd_eu_intW10
  wstd_auth01W11 ~~ rcov3*wstd_eu_intW11
  wstd_auth01W14 ~~ rcov4*wstd_eu_intW14
  
  
  # Estimate the variance and covariance of the random intercepts. 
  RIstd_auth01W ~~ 0.729*RIstd_auth01W
  RIstd_eu_intW ~~ 0.841*RIstd_eu_intW
  RIstd_auth01W ~~ RIstd_eu_intW
  
  # Set variances of within-components at first wave to 1
  wstd_auth01W7 ~~ 1*wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ 1*wstd_eu_intW7 
  
  # Label residual variances
  wstd_auth01W10 ~~ rvx2*wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ rvy2*wstd_eu_intW10 
  wstd_auth01W11 ~~ rvx3*wstd_auth01W11
  wstd_eu_intW11 ~~ rvy3*wstd_eu_intW11 
  wstd_auth01W14 ~~ rvx4*wstd_auth01W14 
  wstd_eu_intW14 ~~ rvy4*wstd_eu_intW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W10 ~~ 0.0058*std_auth01W10
std_auth01W11 ~~ 0.002*std_auth01W11
std_auth01W14 ~~ 0.00725*std_auth01W14
std_eu_intW10 ~~ 0.0028*std_eu_intW10
std_eu_intW11 ~~ 0.00194*std_eu_intW11
std_eu_intW14 ~~ 0*std_eu_intW14

  # Compute correlations of within-components at each wave
  cor2 := a*c + b*d + a*d*cor1 + b*c*cor1 + rcov2
  cor3 := a*c + b*d + a*d*cor2 + b*c*cor2 + rcov3
  cor4 := a*c + b*d + a*d*cor3 + b*c*cor3 + rcov4
 
  # Contrain residual variances of within-components such that variance of each 
  # within-component equals 1
  rvx2 == 1 - (a*a + b*b + 2*a*b*cor1)
  rvy2 == 1 - (c*c + d*d + 2*c*d*cor1)
  rvx3 == 1 - (a*a + b*b + 2*a*b*cor2)
  rvy3 == 1 - (c*c + d*d + 2*c*d*cor2)
  rvx4 == 1 - (a*a + b*b + 2*a*b*cor3)
  rvy4 == 1 - (c*c + d*d + 2*c*d*cor3)
  
  # incorporate tvc - age
wstd_auth01W7 + wstd_eu_intW7 ~ ageW7
wstd_auth01W10 + wstd_eu_intW10 ~ ageW10
wstd_auth01W11 + wstd_eu_intW11 ~ ageW11
wstd_auth01W14 + wstd_eu_intW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc – income
wstd_auth01W7 + wstd_eu_intW7 ~ incomerec7
wstd_auth01W10 + wstd_eu_intW10 ~ incomerec10
wstd_auth01W11 + wstd_eu_intW11 ~ incomerec11
wstd_auth01W14 + wstd_eu_intW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14

# incorporate tvc – social grade
wstd_auth01W7 + wstd_eu_intW7 ~ p_socgradeW7
wstd_auth01W10 + wstd_eu_intW10 ~ p_socgradeW10
wstd_auth01W11 + wstd_eu_intW11 ~ p_socgradeW11
wstd_auth01W14 + wstd_eu_intW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14
'

std_eu_EQ.fit <- lavaan(std_eu_EQ, 
                               data = dat, 
                               estimator = "MLR",
                               missing = 'ML', 
                               meanstructure = T, 
                               int.ov.free = T) 
summary(std_eu_EQ.fit, standardized = T)

#### 7. Model fit - RI-CLPMs -----------------

## Comparison of time-varying and constrained models: 

# Time-varying 
fitmeasures(std_immig_pref_2_tvc.fit)
fitmeasures(std_eu_pref_tvc.fit)

#Time-invariant 
fitmeasures(std_immig_pref_EQ.fit)
fitmeasures(std_eu_EQ.fit)

## Immigration
# Time-varying: CFI = 0.995; AIC = 435,701.424; RMSEA = 0.023; SRMR = 0.06 
# Constrained: CFI = 0.994; AIC = 435,715.118; RMSEA = 0.022; SRMR = 0.061

## EU
# Time-varying: CFI = 0.994; AIC = 440,728.583; RMSEA = 0.023; SRMR = 0.057
# Constrained: 

#### 8. Results - RI-CLPMs ####

standardizedsolution(std_immig_pref_2_tvc.fit)
parameterestimates(std_immig_pref_EQ.fit)

standardizedsolution(std_eu_pref_tvc.fit)
parameterestimates(std_eu_EQ.fit)


#### 9. CLPMs ------

# Immigration 

std_immig_pref_CLPM_tvc <- '
  # Create between components (random intercepts)
  RIauth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wauth01W7 =~ 1*std_auth01W7
  wauth01W10 =~ 1*std_auth01W10
  wauth01W11 =~ 1*std_auth01W11 
  wauth01W14 =~ 1*std_auth01W14

  wim_prefW7 =~ 1*std_im_pref_2W7
  wim_prefW10 =~ 1*std_im_pref_2W10
  wim_prefW11 =~ 1*std_im_pref_2W11
  wim_prefW14 =~ 1*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wauth01W10 + wim_prefW10 ~ wauth01W7 + wim_prefW7
  wauth01W11 + wim_prefW11 ~ wauth01W10 + wim_prefW10
  wauth01W14 + wim_prefW14 ~ wauth01W11 + wim_prefW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wauth01W7 ~~ wim_prefW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wauth01W10 ~~ wim_prefW10
  wauth01W11 ~~ wim_prefW11
  wauth01W14 ~~ wim_prefW14
  
  
  # Fix the variance and covariance of the random intercepts to zero 
  RIauth01W ~~ 0*RIauth01W
  RIim_prefW ~~ 0*RIim_prefW
  RIauth01W ~~ 0*RIim_prefW

  # Estimate the (residual) variance of the within-person centered variables.
  wauth01W7 ~~ wauth01W7 # Variances
  wim_prefW7 ~~ wim_prefW7 
  wauth01W10 ~~ wauth01W10 # Residual variances
  wim_prefW10 ~~ wim_prefW10 
  wauth01W11 ~~ wauth01W11
  wim_prefW11 ~~ wim_prefW11 
  wauth01W14 ~~ wauth01W14 
  wim_prefW14 ~~ wim_prefW14 
  
# Fix error variance of indicators to (1 - alpha)*var
std_auth01W7 ~~ 0.02*std_auth01W7
std_auth01W10 ~~ 0.016*std_auth01W10
std_auth01W11 ~~ 0.018*std_auth01W11
std_auth01W14 ~~ 0.019*std_auth01W14
std_im_pref_2W7 ~~ 0.0147*std_im_pref_2W7
std_im_pref_2W10 ~~ 0.000155*std_im_pref_2W10
std_im_pref_2W11 ~~ 0.000147*std_im_pref_2W11
std_im_pref_2W14 ~~ 0.000512*std_im_pref_2W14

# incorporate tvc - age
wauth01W7 + wim_prefW7 ~ ageW7
wauth01W10 + wim_prefW10 ~ ageW10
wauth01W11 + wim_prefW11 ~ ageW11
wauth01W14 + wim_prefW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc - ethnicity
wauth01W7 + wim_prefW7 ~ ethnic7
wauth01W10 + wim_prefW10 ~ ethnic10
wauth01W11 + wim_prefW11 ~ ethnic11
wauth01W14 + wim_prefW14 ~ ethnic14

ethnic7 ~~ ethnic7
ethnic10 ~~ ethnic10
ethnic11 ~~ ethnic11
ethnic14 ~~ ethnic14

ethnic7 ~~ ethnic10 + ethnic11 + ethnic14
ethnic10 ~~ ethnic11 + ethnic14
ethnic11 ~~ ethnic14

# incorporate tvc – social grade
wauth01W7 + wim_prefW7 ~ p_socgradeW7
wauth01W10 + wim_prefW10 ~ p_socgradeW10
wauth01W11 + wim_prefW11 ~ p_socgradeW11
wauth01W14 + wim_prefW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – education
wauth01W7 + wim_prefW7 ~ higheduc7
wauth01W10 + wim_prefW10 ~ higheduc10
wauth01W11 + wim_prefW11 ~ higheduc11
wauth01W14 + wim_prefW14 ~ higheduc14

higheduc7 ~~ higheduc7
higheduc10 ~~ higheduc10
higheduc11 ~~ higheduc11
higheduc14 ~~ higheduc14

higheduc7 ~~ higheduc10 + higheduc11 + higheduc14
higheduc10 ~~ higheduc11 + higheduc14
higheduc11 ~~ higheduc14

# incorporate tvc – income
wauth01W7 + wim_prefW7 ~ incomerec7
wauth01W10 + wim_prefW10 ~ incomerec10
wauth01W11 + wim_prefW11 ~ incomerec11
wauth01W14 + wim_prefW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14


  '
std_immig_pref_CLPM_tvc.fit <- lavaan(std_immig_pref_CLPM_tvc,
                                      data = dat, 
                                      estimator = "MLR",
                                      missing = 'ML', 
                                      meanstructure = T, 
                                      int.ov.free = T) 
summary(std_immig_pref_CLPM_tvc.fit, standardized = T)

# EU integration 

std_eu_pref_CLPM_tvc <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wstd_eu_intW7 =~ 1*std_eu_intW7
  wstd_eu_intW10 =~ 1*std_eu_intW10
  wstd_eu_intW11 =~ 1*std_eu_intW11
  wstd_eu_intW14 =~ 1*std_eu_intW14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 + wstd_eu_intW10 ~ wstd_auth01W7 + wstd_eu_intW7
  wstd_auth01W11 + wstd_eu_intW11 ~ wstd_auth01W10 + wstd_eu_intW10
  wstd_auth01W14 + wstd_eu_intW14 ~ wstd_auth01W11 + wstd_eu_intW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wstd_eu_intW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wstd_eu_intW10
  wstd_auth01W11 ~~ wstd_eu_intW11
  wstd_auth01W14 ~~ wstd_eu_intW14
  
  
  # Estimate the variance and covariance of the random intercepts. 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIstd_eu_intW ~~ 0*RIstd_eu_intW
  RIstd_auth01W ~~ 0*RIstd_eu_intW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ wstd_eu_intW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ wstd_eu_intW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wstd_eu_intW11 ~~ wstd_eu_intW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wstd_eu_intW14 ~~ wstd_eu_intW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W10 ~~ 0.0058*std_auth01W10
std_auth01W11 ~~ 0.002*std_auth01W11
std_auth01W14 ~~ 0.00725*std_auth01W14
std_eu_intW10 ~~ 0.0028*std_eu_intW10
std_eu_intW11 ~~ 0.00194*std_eu_intW11
std_eu_intW14 ~~ 0*std_eu_intW14

# incorporate tvc - age
wstd_auth01W7 + wstd_eu_intW7 ~ ageW7
wstd_auth01W10 + wstd_eu_intW10 ~ ageW10
wstd_auth01W11 + wstd_eu_intW11 ~ ageW11
wstd_auth01W14 + wstd_eu_intW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc - ethnicity
wstd_auth01W7 + wstd_eu_intW7 ~ ethnic7
wstd_auth01W10 + wstd_eu_intW10 ~ ethnic10
wstd_auth01W11 + wstd_eu_intW11 ~ ethnic11
wstd_auth01W14 + wstd_eu_intW14 ~ ethnic14

ethnic7 ~~ ethnic7
ethnic10 ~~ ethnic10
ethnic11 ~~ ethnic11
ethnic14 ~~ ethnic14

ethnic7 ~~ ethnic10 + ethnic11 + ethnic14
ethnic10 ~~ ethnic11 + ethnic14
ethnic11 ~~ ethnic14

# incorporate tvc – social grade
wstd_auth01W7 + wstd_eu_intW7 ~ p_socgradeW7
wstd_auth01W10 + wstd_eu_intW10 ~ p_socgradeW10
wstd_auth01W11 + wstd_eu_intW11 ~ p_socgradeW11
wstd_auth01W14 + wstd_eu_intW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – education
wstd_auth01W7 + wstd_eu_intW7 ~ higheduc7
wstd_auth01W10 + wstd_eu_intW10 ~ higheduc10
wstd_auth01W11 + wstd_eu_intW11 ~ higheduc11
wstd_auth01W14 + wstd_eu_intW14 ~ higheduc14

higheduc7 ~~ higheduc7
higheduc10 ~~ higheduc10
higheduc11 ~~ higheduc11
higheduc14 ~~ higheduc14

higheduc7 ~~ higheduc10 + higheduc11 + higheduc14
higheduc10 ~~ higheduc11 + higheduc14
higheduc11 ~~ higheduc14

# incorporate tvc – income
wstd_auth01W7 + wstd_eu_intW7 ~ incomerec7
wstd_auth01W10 + wstd_eu_intW10 ~ incomerec10
wstd_auth01W11 + wstd_eu_intW11 ~ incomerec11
wstd_auth01W14 + wstd_eu_intW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14

'

std_eu_pref_CLPM_tvc.fit <- lavaan(std_eu_pref_CLPM_tvc, data = dat, 
                                   estimator = "MLR",
                                   missing = 'ML', meanstructure = T, int.ov.free = T) 
summary(std_eu_pref_CLPM_tvc.fit, standardized = T)


#### 10. CLPMs - equality constraints --------

# Immigration

immig_pref_CLPM_EQ <- '
  # Create between components (random intercepts)
  RIauth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wauth01W7 =~ NA*std_auth01W7
  wauth01W10 =~ NA*std_auth01W10
  wauth01W11 =~ NA*std_auth01W11 
  wauth01W14 =~ NA*std_auth01W14

  wim_prefW7 =~ NA*std_im_pref_2W7
  wim_prefW10 =~ NA*std_im_pref_2W10
  wim_prefW11 =~ NA*std_im_pref_2W11
  wim_prefW14 =~ NA*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wauth01W10 ~ a*wauth01W7 + b*wim_prefW7
  wim_prefW10 ~ c*wauth01W7 + d*wim_prefW7
  wauth01W11 ~ a*wauth01W10 + b*wim_prefW10
  wim_prefW11 ~ c*wauth01W10 + d*wim_prefW10
  wauth01W14 ~ a*wauth01W11 + b*wim_prefW11
  wim_prefW14 ~ c*wauth01W11 + d*wim_prefW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wauth01W7 ~~ cor1*wim_prefW7 # Covariance
  
  # Label residual covariances
  wauth01W10 ~~ rcov2*wim_prefW10
  wauth01W11 ~~ rcov3*wim_prefW11
  wauth01W14 ~~ rcov4*wim_prefW14
  
  
  # Fix the variance and covariance of the random intercepts to zero 
  RIauth01W ~~ 0*RIauth01W
  RIim_prefW ~~ 0*RIim_prefW
  RIauth01W ~~ 0*RIim_prefW

  # Set variances of within-components at first wave to 1
  wauth01W7 ~~ 1*wauth01W7 # Variances
  wim_prefW7 ~~ 1*wim_prefW7
  
  wauth01W10 ~~ rvx2*wauth01W10 # Residual variances
  wim_prefW10 ~~ rvy2*wim_prefW10 
  wauth01W11 ~~ rvx3*wauth01W11
  wim_prefW11 ~~ rvy3*wim_prefW11 
  wauth01W14 ~~ rvx4*wauth01W14 
  wim_prefW14 ~~ rvy4*wim_prefW14 
  
# Fix error variance of indicators to (1 - alpha)*var
std_auth01W7 ~~ 0.02*std_auth01W7
std_auth01W10 ~~ 0.016*std_auth01W10
std_auth01W11 ~~ 0.018*std_auth01W11
std_auth01W14 ~~ 0.019*std_auth01W14
std_im_pref_2W7 ~~ 0.0147*std_im_pref_2W7
std_im_pref_2W10 ~~ 0.000155*std_im_pref_2W10
std_im_pref_2W11 ~~ 0.000147*std_im_pref_2W11
std_im_pref_2W14 ~~ 0.000512*std_im_pref_2W14

 # Compute correlations of within-components at each wave
  cor2 := a*c + b*d + a*d*cor1 + b*c*cor1 + rcov2
  cor3 := a*c + b*d + a*d*cor2 + b*c*cor2 + rcov3
  cor4 := a*c + b*d + a*d*cor3 + b*c*cor3 + rcov4
 
  # Contrain residual variances of within-components such that variance of each 
  # within-component equals 1
  rvx2 == 1 - (a*a + b*b + 2*a*b*cor1)
  rvy2 == 1 - (c*c + d*d + 2*c*d*cor1)
  rvx3 == 1 - (a*a + b*b + 2*a*b*cor2)
  rvy3 == 1 - (c*c + d*d + 2*c*d*cor2)
  rvx4 == 1 - (a*a + b*b + 2*a*b*cor3)
  rvy4 == 1 - (c*c + d*d + 2*c*d*cor3)
  
  # incorporate tvc - age
wauth01W7 + wim_prefW7 ~ ageW7
wauth01W10 + wim_prefW10 ~ ageW10
wauth01W11 + wim_prefW11 ~ ageW11
wauth01W14 + wim_prefW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc - ethnicity
wauth01W7 + wim_prefW7 ~ ethnic7
wauth01W10 + wim_prefW10 ~ ethnic10
wauth01W11 + wim_prefW11 ~ ethnic11
wauth01W14 + wim_prefW14 ~ ethnic14

ethnic7 ~~ ethnic7
ethnic10 ~~ ethnic10
ethnic11 ~~ ethnic11
ethnic14 ~~ ethnic14

ethnic7 ~~ ethnic10 + ethnic11 + ethnic14
ethnic10 ~~ ethnic11 + ethnic14
ethnic11 ~~ ethnic14

# incorporate tvc – social grade
wauth01W7 + wim_prefW7 ~ p_socgradeW7
wauth01W10 + wim_prefW10 ~ p_socgradeW10
wauth01W11 + wim_prefW11 ~ p_socgradeW11
wauth01W14 + wim_prefW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – education
wauth01W7 + wim_prefW7 ~ higheduc7
wauth01W10 + wim_prefW10 ~ higheduc10
wauth01W11 + wim_prefW11 ~ higheduc11
wauth01W14 + wim_prefW14 ~ higheduc14

higheduc7 ~~ higheduc7
higheduc10 ~~ higheduc10
higheduc11 ~~ higheduc11
higheduc14 ~~ higheduc14

higheduc7 ~~ higheduc10 + higheduc11 + higheduc14
higheduc10 ~~ higheduc11 + higheduc14
higheduc11 ~~ higheduc14

# incorporate tvc – income
wauth01W7 + wim_prefW7 ~ incomerec7
wauth01W10 + wim_prefW10 ~ incomerec10
wauth01W11 + wim_prefW11 ~ incomerec11
wauth01W14 + wim_prefW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14

  '
immig_pref_CLPM_EQ.fit <- lavaan(immig_pref_CLPM_EQ,
                                     data = dat, 
                                     estimator = "MLR",
                                     missing = 'ML', 
                                     meanstructure = T, 
                                     int.ov.free = T) 
summary(immig_pref_CLPM_EQ.fit, standardized = T)

# EU integration


eu_CLPM_EQ <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ NA*std_auth01W7
  wstd_auth01W10 =~ NA*std_auth01W10
  wstd_auth01W11 =~ NA*std_auth01W11 
  wstd_auth01W14 =~ NA*std_auth01W14

  wstd_eu_intW7 =~ NA*std_eu_intW7
  wstd_eu_intW10 =~ NA*std_eu_intW10
  wstd_eu_intW11 =~ NA*std_eu_intW11
  wstd_eu_intW14 =~ NA*std_eu_intW14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 ~ a*wstd_auth01W7 + b*wstd_eu_intW7
  wstd_eu_intW10 ~ c*wstd_auth01W7 + d*wstd_eu_intW7
  wstd_auth01W11 ~ a*wstd_auth01W10 + b*wstd_eu_intW10
  wstd_eu_intW11 ~ c*wstd_auth01W10 + d*wstd_eu_intW10
  wstd_auth01W14 ~ a*wstd_auth01W11 + b*wstd_eu_intW11
  wstd_eu_intW14 ~ c*wstd_auth01W11 + d*wstd_eu_intW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ cor1*wstd_eu_intW7 # Covariance
  
  # Label residual covariances
  wstd_auth01W10 ~~ rcov2*wstd_eu_intW10
  wstd_auth01W11 ~~ rcov3*wstd_eu_intW11
  wstd_auth01W14 ~~ rcov4*wstd_eu_intW14
  
  
  # Fix the variance and covariance of the random intercepts to zero
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIstd_eu_intW ~~ 0*RIstd_eu_intW
  RIstd_auth01W ~~ 0*RIstd_eu_intW
  
  # Set variances of within-components at first wave to 1
  wstd_auth01W7 ~~ 1*wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ 1*wstd_eu_intW7 
  
  # Label residual variances
  wstd_auth01W10 ~~ rvx2*wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ rvy2*wstd_eu_intW10 
  wstd_auth01W11 ~~ rvx3*wstd_auth01W11
  wstd_eu_intW11 ~~ rvy3*wstd_eu_intW11 
  wstd_auth01W14 ~~ rvx4*wstd_auth01W14 
  wstd_eu_intW14 ~~ rvy4*wstd_eu_intW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W10 ~~ 0.0058*std_auth01W10
std_auth01W11 ~~ 0.002*std_auth01W11
std_auth01W14 ~~ 0.00725*std_auth01W14
std_eu_intW10 ~~ 0.0028*std_eu_intW10
std_eu_intW11 ~~ 0.00194*std_eu_intW11
std_eu_intW14 ~~ 0*std_eu_intW14

  # Compute correlations of within-components at each wave
  cor2 := a*c + b*d + a*d*cor1 + b*c*cor1 + rcov2
  cor3 := a*c + b*d + a*d*cor2 + b*c*cor2 + rcov3
  cor4 := a*c + b*d + a*d*cor3 + b*c*cor3 + rcov4
 
  # Contrain residual variances of within-components such that variance of each 
  # within-component equals 1
  rvx2 == 1 - (a*a + b*b + 2*a*b*cor1)
  rvy2 == 1 - (c*c + d*d + 2*c*d*cor1)
  rvx3 == 1 - (a*a + b*b + 2*a*b*cor2)
  rvy3 == 1 - (c*c + d*d + 2*c*d*cor2)
  rvx4 == 1 - (a*a + b*b + 2*a*b*cor3)
  rvy4 == 1 - (c*c + d*d + 2*c*d*cor3)

# incorporate tvc - age
wstd_auth01W7 + wstd_eu_intW7 ~ ageW7
wstd_auth01W10 + wstd_eu_intW10 ~ ageW10
wstd_auth01W11 + wstd_eu_intW11 ~ ageW11
wstd_auth01W14 + wstd_eu_intW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc - ethnicity
wstd_auth01W7 + wstd_eu_intW7 ~ ethnic7
wstd_auth01W10 + wstd_eu_intW10 ~ ethnic10
wstd_auth01W11 + wstd_eu_intW11 ~ ethnic11
wstd_auth01W14 + wstd_eu_intW14 ~ ethnic14

ethnic7 ~~ ethnic7
ethnic10 ~~ ethnic10
ethnic11 ~~ ethnic11
ethnic14 ~~ ethnic14

ethnic7 ~~ ethnic10 + ethnic11 + ethnic14
ethnic10 ~~ ethnic11 + ethnic14
ethnic11 ~~ ethnic14

# incorporate tvc – social grade
wstd_auth01W7 + wstd_eu_intW7 ~ p_socgradeW7
wstd_auth01W10 + wstd_eu_intW10 ~ p_socgradeW10
wstd_auth01W11 + wstd_eu_intW11 ~ p_socgradeW11
wstd_auth01W14 + wstd_eu_intW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – education
wstd_auth01W7 + wstd_eu_intW7 ~ higheduc7
wstd_auth01W10 + wstd_eu_intW10 ~ higheduc10
wstd_auth01W11 + wstd_eu_intW11 ~ higheduc11
wstd_auth01W14 + wstd_eu_intW14 ~ higheduc14

higheduc7 ~~ higheduc7
higheduc10 ~~ higheduc10
higheduc11 ~~ higheduc11
higheduc14 ~~ higheduc14

higheduc7 ~~ higheduc10 + higheduc11 + higheduc14
higheduc10 ~~ higheduc11 + higheduc14
higheduc11 ~~ higheduc14

# incorporate tvc – income
wstd_auth01W7 + wstd_eu_intW7 ~ incomerec7
wstd_auth01W10 + wstd_eu_intW10 ~ incomerec10
wstd_auth01W11 + wstd_eu_intW11 ~ incomerec11
wstd_auth01W14 + wstd_eu_intW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14

'

eu_CLPM_EQ.fit <- lavaan(eu_CLPM_EQ, 
                             data = dat, 
                             estimator = "MLR",
                             missing = 'ML', 
                             meanstructure = T, 
                             int.ov.free = T) 
summary(eu_CLPM_EQ.fit, standardized = T)


#### 11. Model fit - CLPMs ----------


# Comparison of time-varying and constrained models: 

# Time-varying 
fitmeasures(std_immig_pref_CLPM_tvc.fit)
fitmeasures(std_eu_pref_CLPM_tvc.fit)

#Time-invariant 
fitmeasures(immig_pref_CLPM_EQ.fit)
fitmeasures(eu_CLPM_EQ.fit)

## Immigration
# Time-varying: CFI = 0.98; AIC = 411812.368; RMSEA = 0.035; SRMR = 0.09 
# Constrained: CFI = 0.98; AIC =411877.677 ; RMSEA = 0.035; SRMR = 0.09

## EU
# Time-varying: CFI = 0.979; AIC = 417076.560; RMSEA = 0.035; SRMR = 0.09
# Constrained: CFI = 0.978; AIC = 417294.930; RMSEA = 0.035; SRMR = 0.09

#### 12. Results - CLPMs ----------

standardizedsolution(std_immig_pref_CLPM_tvc.fit)
parameterestimates(immig_pref_CLPM_EQ.fit)

standardizedsolution(std_eu_pref_CLPM_tvc.fit)
parameterestimates(eu_CLPM_EQ.fit)

#### 13. Multiple-group RI-CLPMs - immigration -------------

# New data frame for moderation analysis

head(dat)

att_dat <- subset(dat, polatt_group == 0 | polatt_group == 1)

# First fit a model in which differences between the groups are left
# unconstrained 

std_immig_attU <- '
 # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wim_prefW7 =~ 1*std_im_pref_2W7
  wim_prefW10 =~ 1*std_im_pref_2W10
  wim_prefW11 =~ 1*std_im_pref_2W11
  wim_prefW14 =~ 1*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 + wim_prefW10 ~ wstd_auth01W7 + wim_prefW7
  wstd_auth01W11 + wim_prefW11 ~ wstd_auth01W10 + wim_prefW10
  wstd_auth01W14 + wim_prefW14 ~ wstd_auth01W11 + wim_prefW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wim_prefW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wim_prefW10
  wstd_auth01W11 ~~ wim_prefW11
  wstd_auth01W14 ~~ wim_prefW14
  
  
  # Estimate the variance and covariance of the random intercepts. 
  RIstd_auth01W ~~ RIstd_auth01W
  RIim_prefW ~~ RIim_prefW
  RIstd_auth01W ~~ RIim_prefW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wim_prefW7 ~~ wim_prefW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wim_prefW10 ~~ wim_prefW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wim_prefW11 ~~ wim_prefW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wim_prefW14 ~~ wim_prefW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W7 ~~ 0.02*std_auth01W7
std_auth01W10 ~~ 0.016*std_auth01W10
std_auth01W11 ~~ 0.018*std_auth01W11
std_auth01W14 ~~ 0.019*std_auth01W14
std_im_pref_2W7 ~~ 0.0147*std_im_pref_2W7
std_im_pref_2W10 ~~ 0.000155*std_im_pref_2W10
std_im_pref_2W11 ~~ 0.000147*std_im_pref_2W11
std_im_pref_2W14 ~~ 0.000512*std_im_pref_2W14

  '
std_immig_attU.fit <- lavaan(std_immig_attU, 
                             data = att_dat, 
                             group = "polatt_group",
                             estimator = "MLR",
                             missing = 'ML',
                             meanstructure = T, 
                             int.ov.free = T) 
summary(std_immig_attU.fit, standardized = T)

# Then fit a model in which the two groups are constrained to equality

std_immig_attC <- '
 # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wim_prefW7 =~ 1*std_im_pref_2W7
  wim_prefW10 =~ 1*std_im_pref_2W10
  wim_prefW11 =~ 1*std_im_pref_2W11
  wim_prefW14 =~ 1*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 ~ c(a1, a1)*wstd_auth01W7 + c(b1, b1)*wim_prefW7
  wim_prefW10 ~ c(c1, c1)*wstd_auth01W7 + c(d1, d1)*wim_prefW7
  wstd_auth01W11 ~ c(a2, a2)*wstd_auth01W10 + c(b2, b2)*wim_prefW10
  wim_prefW11 ~ c(c2, c2)*wstd_auth01W10 + c(d2, d2)*wim_prefW10
  wstd_auth01W14 ~ c(a3, a3)*wstd_auth01W11 + c(b3, b3)*wim_prefW11
  wim_prefW14 ~ c(c3, c3)*wstd_auth01W11 + c(d3, d3)*wim_prefW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wim_prefW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wim_prefW10
  wstd_auth01W11 ~~ wim_prefW11
  wstd_auth01W14 ~~ wim_prefW14
  
  
  # Estimate the variance and covariance of the random intercepts. 
  RIstd_auth01W ~~ RIstd_auth01W
  RIim_prefW ~~ RIim_prefW
  RIstd_auth01W ~~ RIim_prefW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wim_prefW7 ~~ wim_prefW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wim_prefW10 ~~ wim_prefW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wim_prefW11 ~~ wim_prefW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wim_prefW14 ~~ wim_prefW14 

  '
std_immig_attC.fit <- lavaan(std_immig_attC, 
                             data = att_dat, 
                             group = "polatt_group",
                             estimator = "MLR",
                             missing = 'ML',
                             meanstructure = T, 
                             int.ov.free = T) 
summary(std_immig_attC.fit, standardized = T)

# Then compare model fit using an ANOVA

anova(std_immig_attU.fit, std_immig_attC.fit)

# p > 0.05, indicating that the unconstrained model does not fit the data 
# better and thus that there are no differences between the groups

#### 14. Multiple group RI-CLPMs - EU -----------

# Begin with an unconstrained model 

std_eu_attU <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wstd_eu_intW7 =~ 1*std_eu_intW7
  wstd_eu_intW10 =~ 1*std_eu_intW10
  wstd_eu_intW11 =~ 1*std_eu_intW11
  wstd_eu_intW14 =~ 1*std_eu_intW14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 + wstd_eu_intW10 ~ wstd_auth01W7 + wstd_eu_intW7
  wstd_auth01W11 + wstd_eu_intW11 ~ wstd_auth01W10 + wstd_eu_intW10
  wstd_auth01W14 + wstd_eu_intW14 ~ wstd_auth01W11 + wstd_eu_intW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wstd_eu_intW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wstd_eu_intW10
  wstd_auth01W11 ~~ wstd_eu_intW11
  wstd_auth01W14 ~~ wstd_eu_intW14
  
  
  # Estimate the variance and covariance of the random intercepts. 
  RIstd_auth01W ~~ RIstd_auth01W
  RIstd_eu_intW ~~ RIstd_eu_intW
  RIstd_auth01W ~~ RIstd_eu_intW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ wstd_eu_intW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ wstd_eu_intW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wstd_eu_intW11 ~~ wstd_eu_intW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wstd_eu_intW14 ~~ wstd_eu_intW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W10 ~~ 0.0058*std_auth01W10
std_auth01W11 ~~ 0.002*std_auth01W11
std_auth01W14 ~~ 0.00725*std_auth01W14
std_eu_intW10 ~~ 0.0028*std_eu_intW10
std_eu_intW11 ~~ 0.00194*std_eu_intW11
std_eu_intW14 ~~ 0*std_eu_intW14


'

std_eu_attU.fit <- lavaan(std_eu_attU, 
                          estimator = "MLR",
                          data = att_dat,
                          group = "polatt_group",
                          missing = 'ML', 
                          meanstructure = T, 
                          int.ov.free = T) 
summary(std_eu_attU.fit, standardized = T)

# Then fit a constrained model 

std_eu_attC <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wstd_eu_intW7 =~ 1*std_eu_intW7
  wstd_eu_intW10 =~ 1*std_eu_intW10
  wstd_eu_intW11 =~ 1*std_eu_intW11
  wstd_eu_intW14 =~ 1*std_eu_intW14
  
  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 ~ c(a1, a1)*wstd_auth01W7 + c(b1, b1)*wstd_eu_intW7
  wstd_eu_intW10 ~ c(c1, c1)*wstd_auth01W7 + c(d1, d1)*wstd_eu_intW7
  wstd_auth01W11 ~ c(a2, a2)*wstd_auth01W10 + c(b2, b2)*wstd_eu_intW10
  wstd_eu_intW11 ~ c(c2, c2)*wstd_auth01W10 + c(d2, d2)*wstd_eu_intW10
  wstd_auth01W14 ~ c(a3, a3)*wstd_auth01W11 + c(b3, b3)*wstd_eu_intW11
  wstd_eu_intW14 ~ c(c3, c3)*wstd_auth01W11 + c(d3, d3)*wstd_eu_intW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wstd_eu_intW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wstd_eu_intW10
  wstd_auth01W11 ~~ wstd_eu_intW11
  wstd_auth01W14 ~~ wstd_eu_intW14
  
  
  # Estimate the variance and covariance of the random intercepts. 
  RIstd_auth01W ~~ RIstd_auth01W
  RIstd_eu_intW ~~ RIstd_eu_intW
  RIstd_auth01W ~~ RIstd_eu_intW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ wstd_eu_intW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ wstd_eu_intW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wstd_eu_intW11 ~~ wstd_eu_intW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wstd_eu_intW14 ~~ wstd_eu_intW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W10 ~~ 0.0058*std_auth01W10
std_auth01W11 ~~ 0.002*std_auth01W11
std_auth01W14 ~~ 0.00725*std_auth01W14
std_eu_intW10 ~~ 0.0028*std_eu_intW10
std_eu_intW11 ~~ 0.00194*std_eu_intW11
std_eu_intW14 ~~ 0*std_eu_intW14
'

std_eu_attC.fit <- lavaan(std_eu_attC, 
                          estimator = "MLR",
                          data = att_dat,
                          group = "polatt_group",
                          missing = 'ML', 
                          meanstructure = T, 
                          int.ov.free = T) 
summary(std_eu_attC.fit, standardized = T)

# ANOVA for model fit comparison

anova(std_eu_attC.fit, std_eu_attU.fit)

# Again, p > 0.05

#### 15. Multiple-group CLPMs - immigration -------------

std_immig_attU_CLPM <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wim_prefW7 =~ 1*std_im_pref_2W7
  wim_prefW10 =~ 1*std_im_pref_2W10
  wim_prefW11 =~ 1*std_im_pref_2W11
  wim_prefW14 =~ 1*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 + wim_prefW10 ~ wstd_auth01W7 + wim_prefW7
  wstd_auth01W11 + wim_prefW11 ~ wstd_auth01W10 + wim_prefW10
  wstd_auth01W14 + wim_prefW14 ~ wstd_auth01W11 + wim_prefW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wim_prefW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wim_prefW10
  wstd_auth01W11 ~~ wim_prefW11
  wstd_auth01W14 ~~ wim_prefW14
  
  
  # Constrain the variance and covariance of the random intercepts to 0 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIim_prefW ~~ 0*RIim_prefW
  RIstd_auth01W ~~ 0*RIim_prefW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wim_prefW7 ~~ wim_prefW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wim_prefW10 ~~ wim_prefW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wim_prefW11 ~~ wim_prefW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wim_prefW14 ~~ wim_prefW14 
  
# Fix error variance of indicators to (1 - alpha)*var
std_auth01W7 ~~ 0.02*std_auth01W7
std_auth01W10 ~~ 0.016*std_auth01W10
std_auth01W11 ~~ 0.018*std_auth01W11
std_auth01W14 ~~ 0.019*std_auth01W14
std_im_pref_2W7 ~~ 0.0147*std_im_pref_2W7
std_im_pref_2W10 ~~ 0.000155*std_im_pref_2W10
std_im_pref_2W11 ~~ 0.000147*std_im_pref_2W11
std_im_pref_2W14 ~~ 0.000512*std_im_pref_2W14
  '
std_immig_attU_CLPM.fit <- lavaan(std_immig_attU_CLPM, 
                                  data = att_dat, 
                                  group = "polatt_group",
                                  estimator = "MLR",
                                  missing = 'ML',
                                  meanstructure = T, 
                                  int.ov.free = T) 
summary(std_immig_attU_CLPM.fit, standardized = T)

std_immig_attC_CLPM <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wim_prefW7 =~ 1*std_im_pref_2W7
  wim_prefW10 =~ 1*std_im_pref_2W10
  wim_prefW11 =~ 1*std_im_pref_2W11
  wim_prefW14 =~ 1*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 ~ c(a1, a1)*wstd_auth01W7 + c(b1, b1)*wim_prefW7
  wim_prefW10 ~ c(c1, c1)*wstd_auth01W7 + c(d1, d1)*wim_prefW7
  wstd_auth01W11 ~ c(a2, a2)*wstd_auth01W10 + c(b2, b2)*wim_prefW10
  wim_prefW11 ~ c(c2, c2)*wstd_auth01W10 + c(d2, d2)*wim_prefW10
  wstd_auth01W14 ~ c(a3, a3)*wstd_auth01W11 + c(b3, b3)*wim_prefW11
  wim_prefW14 ~ c(c3, c3)*wstd_auth01W11 + c(d3, d3)*wim_prefW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wim_prefW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wim_prefW10
  wstd_auth01W11 ~~ wim_prefW11
  wstd_auth01W14 ~~ wim_prefW14
  
  
  # Constrain the variance and covariance of the random intercepts to 0 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIim_prefW ~~ 0*RIim_prefW
  RIstd_auth01W ~~ 0*RIim_prefW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wim_prefW7 ~~ wim_prefW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wim_prefW10 ~~ wim_prefW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wim_prefW11 ~~ wim_prefW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wim_prefW14 ~~ wim_prefW14 
  
# Fix error variance of indicators to (1 - alpha)*var
std_auth01W7 ~~ 0.02*std_auth01W7
std_auth01W10 ~~ 0.016*std_auth01W10
std_auth01W11 ~~ 0.018*std_auth01W11
std_auth01W14 ~~ 0.019*std_auth01W14
std_im_pref_2W7 ~~ 0.0147*std_im_pref_2W7
std_im_pref_2W10 ~~ 0.000155*std_im_pref_2W10
std_im_pref_2W11 ~~ 0.000147*std_im_pref_2W11
std_im_pref_2W14 ~~ 0.000512*std_im_pref_2W14
  '
std_immig_attC_CLPM.fit <- lavaan(std_immig_attC_CLPM, 
                                  data = att_dat, 
                                  group = "polatt_group",
                                  estimator = "MLR",
                                  missing = 'ML',
                                  meanstructure = T, 
                                  int.ov.free = T) 
summary(std_immig_attC_CLPM.fit, standardized = T)

# ANOVA for model fit comparison

anova(std_immig_attU_CLPM.fit, std_immig_attC_CLPM.fit)

# p < 0.05, indicating that the two groups differ. 

## Now, subset the data frame into two data frames - one for the high engagement
# group and one for the low engagement group. This is necessary in order to fit 
# models with equality constraints imposed on the parameters over time. 

att_dat_high <- subset(dat, polatt_group == 1)

dim(att_dat_high)

att_dat_low <- subset(dat, polatt_group == 0)

dim(att_dat_low)

## Highly engaged group 

# with equality constraints 

std_immig_CLPM_H <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ NA*std_auth01W7
  wstd_auth01W10 =~ NA*std_auth01W10
  wstd_auth01W11 =~ NA*std_auth01W11 
  wstd_auth01W14 =~ NA*std_auth01W14

  wim_prefW7 =~ NA*std_im_pref_2W7
  wim_prefW10 =~ NA*std_im_pref_2W10
  wim_prefW11 =~ NA*std_im_pref_2W11
  wim_prefW14 =~ NA*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 ~ a*wstd_auth01W7 + b*wim_prefW7
  wim_prefW10 ~ c*wstd_auth01W7 + d*wim_prefW7
  wstd_auth01W11 ~ a*wstd_auth01W10 + b*wim_prefW10
  wim_prefW11 ~ c*wstd_auth01W10 + d*wim_prefW10
  wstd_auth01W14 ~ a*wstd_auth01W11 + b*wim_prefW11
  wim_prefW14 ~ c*wstd_auth01W11 + d*wim_prefW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ cor1*wim_prefW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ rcov2*wim_prefW10
  wstd_auth01W11 ~~ rcov3*wim_prefW11
  wstd_auth01W14 ~~ rcov4*wim_prefW14
  
  
  # Constrain the variance and covariance of the random intercepts to 0 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIim_prefW ~~ 0*RIim_prefW
  RIstd_auth01W ~~ 0*RIim_prefW

  # Set variances of within-components at first wave to 1
  wstd_auth01W7 ~~ 1*wstd_auth01W7 # Variances
  wim_prefW7 ~~ 1*wim_prefW7 
  
  
  wstd_auth01W10 ~~ rvx2*wstd_auth01W10 # Residual variances
  wim_prefW10 ~~ rvy2*wim_prefW10 
  wstd_auth01W11 ~~ rvx3*wstd_auth01W11
  wim_prefW11 ~~ rvy3*wim_prefW11 
  wstd_auth01W14 ~~ rvx4*wstd_auth01W14 
  wim_prefW14 ~~ rvy4*wim_prefW14 
  
# Fix error variance of indicators to (1 - alpha)*var
std_auth01W7 ~~ 0.02*std_auth01W7
std_auth01W10 ~~ 0.016*std_auth01W10
std_auth01W11 ~~ 0.018*std_auth01W11
std_auth01W14 ~~ 0.019*std_auth01W14
std_im_pref_2W7 ~~ 0.0147*std_im_pref_2W7
std_im_pref_2W10 ~~ 0.000155*std_im_pref_2W10
std_im_pref_2W11 ~~ 0.000147*std_im_pref_2W11
std_im_pref_2W14 ~~ 0.000512*std_im_pref_2W14

 # Compute correlations of within-components at each wave
  cor2 := a*c + b*d + a*d*cor1 + b*c*cor1 + rcov2
  cor3 := a*c + b*d + a*d*cor2 + b*c*cor2 + rcov3
  cor4 := a*c + b*d + a*d*cor3 + b*c*cor3 + rcov4
 
  # Contrain residual variances of within-components such that variance of each 
  # within-component equals 1
  rvx2 == 1 - (a*a + b*b + 2*a*b*cor1)
  rvy2 == 1 - (c*c + d*d + 2*c*d*cor1)
  rvx3 == 1 - (a*a + b*b + 2*a*b*cor2)
  rvy3 == 1 - (c*c + d*d + 2*c*d*cor2)
  rvx4 == 1 - (a*a + b*b + 2*a*b*cor3)
  rvy4 == 1 - (c*c + d*d + 2*c*d*cor3)
  
  # incorporate tvc - age
wstd_auth01W7 + wim_prefW7 ~ ageW7
wstd_auth01W10 + wim_prefW10 ~ ageW10
wstd_auth01W11 + wim_prefW11 ~ ageW11
wstd_auth01W14 + wim_prefW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc - ethnicity
wstd_auth01W7 + wim_prefW7 ~ ethnic7
wstd_auth01W10 + wim_prefW10 ~ ethnic10
wstd_auth01W11 + wim_prefW11 ~ ethnic11
wstd_auth01W14 + wim_prefW14 ~ ethnic14

ethnic7 ~~ ethnic7
ethnic10 ~~ ethnic10
ethnic11 ~~ ethnic11
ethnic14 ~~ ethnic14

ethnic7 ~~ ethnic10 + ethnic11 + ethnic14
ethnic10 ~~ ethnic11 + ethnic14
ethnic11 ~~ ethnic14

# incorporate tvc – social grade
wstd_auth01W7 + wim_prefW7 ~ p_socgradeW7
wstd_auth01W10 + wim_prefW10 ~ p_socgradeW10
wstd_auth01W11 + wim_prefW11 ~ p_socgradeW11
wstd_auth01W14 + wim_prefW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – education
wstd_auth01W7 + wim_prefW7 ~ higheduc7
wstd_auth01W10 + wim_prefW10 ~ higheduc10
wstd_auth01W11 + wim_prefW11 ~ higheduc11
wstd_auth01W14 + wim_prefW14 ~ higheduc14

higheduc7 ~~ higheduc7
higheduc10 ~~ higheduc10
higheduc11 ~~ higheduc11
higheduc14 ~~ higheduc14

higheduc7 ~~ higheduc10 + higheduc11 + higheduc14
higheduc10 ~~ higheduc11 + higheduc14
higheduc11 ~~ higheduc14

# incorporate tvc – income
wstd_auth01W7 + wim_prefW7 ~ incomerec7
wstd_auth01W10 + wim_prefW10 ~ incomerec10
wstd_auth01W11 + wim_prefW11 ~ incomerec11
wstd_auth01W14 + wim_prefW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14


  '
std_immig_CLPM_H.fit <- lavaan(std_immig_CLPM_H, 
                                   data = att_dat_high,
                                   estimator = "MLR",
                                   missing = 'ML',
                                   meanstructure = T, 
                                   int.ov.free = T) 
summary(std_immig_CLPM_H.fit, standardized = T)

# Without equality constraints

std_immig_CLPM_H_noEQ <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wim_prefW7 =~ 1*std_im_pref_2W7
  wim_prefW10 =~ 1*std_im_pref_2W10
  wim_prefW11 =~ 1*std_im_pref_2W11
  wim_prefW14 =~ 1*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 + wim_prefW10 ~ wstd_auth01W7 + wim_prefW7
  wstd_auth01W11 + wim_prefW11 ~ wstd_auth01W10 + wim_prefW10
  wstd_auth01W14 + wim_prefW14 ~ wstd_auth01W11 + wim_prefW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wim_prefW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wim_prefW10
  wstd_auth01W11 ~~ wim_prefW11
  wstd_auth01W14 ~~ wim_prefW14
  
  
  # Constrain the variance and covariance of the random intercepts to 0 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIim_prefW ~~ 0*RIim_prefW
  RIstd_auth01W ~~ 0*RIim_prefW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wim_prefW7 ~~ wim_prefW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wim_prefW10 ~~ wim_prefW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wim_prefW11 ~~ wim_prefW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wim_prefW14 ~~ wim_prefW14 
  
# Fix error variance of indicators to (1 - alpha)*var
std_auth01W7 ~~ 0.02*std_auth01W7
std_auth01W10 ~~ 0.016*std_auth01W10
std_auth01W11 ~~ 0.018*std_auth01W11
std_auth01W14 ~~ 0.019*std_auth01W14
std_im_pref_2W7 ~~ 0.0147*std_im_pref_2W7
std_im_pref_2W10 ~~ 0.000155*std_im_pref_2W10
std_im_pref_2W11 ~~ 0.000147*std_im_pref_2W11
std_im_pref_2W14 ~~ 0.000512*std_im_pref_2W14

  
  # incorporate tvc - age
wstd_auth01W7 + wim_prefW7 ~ ageW7
wstd_auth01W10 + wim_prefW10 ~ ageW10
wstd_auth01W11 + wim_prefW11 ~ ageW11
wstd_auth01W14 + wim_prefW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc - ethnicity
wstd_auth01W7 + wim_prefW7 ~ ethnic7
wstd_auth01W10 + wim_prefW10 ~ ethnic10
wstd_auth01W11 + wim_prefW11 ~ ethnic11
wstd_auth01W14 + wim_prefW14 ~ ethnic14

ethnic7 ~~ ethnic7
ethnic10 ~~ ethnic10
ethnic11 ~~ ethnic11
ethnic14 ~~ ethnic14

ethnic7 ~~ ethnic10 + ethnic11 + ethnic14
ethnic10 ~~ ethnic11 + ethnic14
ethnic11 ~~ ethnic14

# incorporate tvc – social grade
wstd_auth01W7 + wim_prefW7 ~ p_socgradeW7
wstd_auth01W10 + wim_prefW10 ~ p_socgradeW10
wstd_auth01W11 + wim_prefW11 ~ p_socgradeW11
wstd_auth01W14 + wim_prefW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – education
wstd_auth01W7 + wim_prefW7 ~ higheduc7
wstd_auth01W10 + wim_prefW10 ~ higheduc10
wstd_auth01W11 + wim_prefW11 ~ higheduc11
wstd_auth01W14 + wim_prefW14 ~ higheduc14

higheduc7 ~~ higheduc7
higheduc10 ~~ higheduc10
higheduc11 ~~ higheduc11
higheduc14 ~~ higheduc14

higheduc7 ~~ higheduc10 + higheduc11 + higheduc14
higheduc10 ~~ higheduc11 + higheduc14
higheduc11 ~~ higheduc14

# incorporate tvc – income
wstd_auth01W7 + wim_prefW7 ~ incomerec7
wstd_auth01W10 + wim_prefW10 ~ incomerec10
wstd_auth01W11 + wim_prefW11 ~ incomerec11
wstd_auth01W14 + wim_prefW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14
  '
std_immig_CLPM_H_noEQ.fit <- lavaan(std_immig_CLPM_H_noEQ, 
                                        data = att_dat_high,
                                        estimator = "MLR",
                                        missing = 'ML',
                                        meanstructure = T, 
                                        int.ov.free = T) 
summary(std_immig_CLPM_H_noEQ.fit, standardized = T)

## Less engaged group 

# With equality constraints

std_immig_CLPM_L <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ NA*std_auth01W7
  wstd_auth01W10 =~ NA*std_auth01W10
  wstd_auth01W11 =~ NA*std_auth01W11 
  wstd_auth01W14 =~ NA*std_auth01W14

  wim_prefW7 =~ NA*std_im_pref_2W7
  wim_prefW10 =~ NA*std_im_pref_2W10
  wim_prefW11 =~ NA*std_im_pref_2W11
  wim_prefW14 =~ NA*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 ~ a*wstd_auth01W7 + b*wim_prefW7
  wim_prefW10 ~ c*wstd_auth01W7 + d*wim_prefW7
  wstd_auth01W11 ~ a*wstd_auth01W10 + b*wim_prefW10
  wim_prefW11 ~ c*wstd_auth01W10 + d*wim_prefW10
  wstd_auth01W14 ~ a*wstd_auth01W11 + b*wim_prefW11
  wim_prefW14 ~ c*wstd_auth01W11 + d*wim_prefW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ cor1*wim_prefW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ rcov2*wim_prefW10
  wstd_auth01W11 ~~ rcov3*wim_prefW11
  wstd_auth01W14 ~~ rcov4*wim_prefW14
  
  
  # Constrain the variance and covariance of the random intercepts to 0 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIim_prefW ~~ 0*RIim_prefW
  RIstd_auth01W ~~ 0*RIim_prefW

  # Set variances of within-components at first wave to 1
  wstd_auth01W7 ~~ 1*wstd_auth01W7 # Variances
  wim_prefW7 ~~ 1*wim_prefW7 
  
  
  wstd_auth01W10 ~~ rvx2*wstd_auth01W10 # Residual variances
  wim_prefW10 ~~ rvy2*wim_prefW10 
  wstd_auth01W11 ~~ rvx3*wstd_auth01W11
  wim_prefW11 ~~ rvy3*wim_prefW11 
  wstd_auth01W14 ~~ rvx4*wstd_auth01W14 
  wim_prefW14 ~~ rvy4*wim_prefW14 
  
# Fix error variance of indicators to (1 - alpha)*var
std_auth01W7 ~~ 0.02*std_auth01W7
std_auth01W10 ~~ 0.016*std_auth01W10
std_auth01W11 ~~ 0.018*std_auth01W11
std_auth01W14 ~~ 0.019*std_auth01W14
std_im_pref_2W7 ~~ 0.0147*std_im_pref_2W7
std_im_pref_2W10 ~~ 0.000155*std_im_pref_2W10
std_im_pref_2W11 ~~ 0.000147*std_im_pref_2W11
std_im_pref_2W14 ~~ 0.000512*std_im_pref_2W14

 # Compute correlations of within-components at each wave
  cor2 := a*c + b*d + a*d*cor1 + b*c*cor1 + rcov2
  cor3 := a*c + b*d + a*d*cor2 + b*c*cor2 + rcov3
  cor4 := a*c + b*d + a*d*cor3 + b*c*cor3 + rcov4
 
  # Contrain residual variances of within-components such that variance of each 
  # within-component equals 1
  rvx2 == 1 - (a*a + b*b + 2*a*b*cor1)
  rvy2 == 1 - (c*c + d*d + 2*c*d*cor1)
  rvx3 == 1 - (a*a + b*b + 2*a*b*cor2)
  rvy3 == 1 - (c*c + d*d + 2*c*d*cor2)
  rvx4 == 1 - (a*a + b*b + 2*a*b*cor3)
  rvy4 == 1 - (c*c + d*d + 2*c*d*cor3)
  
  # incorporate tvc - age
wstd_auth01W7 + wim_prefW7 ~ ageW7
wstd_auth01W10 + wim_prefW10 ~ ageW10
wstd_auth01W11 + wim_prefW11 ~ ageW11
wstd_auth01W14 + wim_prefW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc - ethnicity
wstd_auth01W7 + wim_prefW7 ~ ethnic7
wstd_auth01W10 + wim_prefW10 ~ ethnic10
wstd_auth01W11 + wim_prefW11 ~ ethnic11
wstd_auth01W14 + wim_prefW14 ~ ethnic14

ethnic7 ~~ ethnic7
ethnic10 ~~ ethnic10
ethnic11 ~~ ethnic11
ethnic14 ~~ ethnic14

ethnic7 ~~ ethnic10 + ethnic11 + ethnic14
ethnic10 ~~ ethnic11 + ethnic14
ethnic11 ~~ ethnic14

# incorporate tvc – social grade
wstd_auth01W7 + wim_prefW7 ~ p_socgradeW7
wstd_auth01W10 + wim_prefW10 ~ p_socgradeW10
wstd_auth01W11 + wim_prefW11 ~ p_socgradeW11
wstd_auth01W14 + wim_prefW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – education
wstd_auth01W7 + wim_prefW7 ~ higheduc7
wstd_auth01W10 + wim_prefW10 ~ higheduc10
wstd_auth01W11 + wim_prefW11 ~ higheduc11
wstd_auth01W14 + wim_prefW14 ~ higheduc14

higheduc7 ~~ higheduc7
higheduc10 ~~ higheduc10
higheduc11 ~~ higheduc11
higheduc14 ~~ higheduc14

higheduc7 ~~ higheduc10 + higheduc11 + higheduc14
higheduc10 ~~ higheduc11 + higheduc14
higheduc11 ~~ higheduc14

# incorporate tvc – income
wstd_auth01W7 + wim_prefW7 ~ incomerec7
wstd_auth01W10 + wim_prefW10 ~ incomerec10
wstd_auth01W11 + wim_prefW11 ~ incomerec11
wstd_auth01W14 + wim_prefW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14


  '
std_immig_CLPM_L.fit <- lavaan(std_immig_CLPM_L, 
                                   data = att_dat_low,
                                   estimator = "MLR",
                                   missing = 'ML',
                                   meanstructure = T, 
                                   int.ov.free = T) 
summary(std_immig_CLPM_L.fit, standardized = T)

# Without equality constraints

std_immig_CLPM_L_noEQ <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIim_prefW =~ 1*std_im_pref_2W7 + 1*std_im_pref_2W10 + 1*std_im_pref_2W11 + 1*std_im_pref_2W14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wim_prefW7 =~ 1*std_im_pref_2W7
  wim_prefW10 =~ 1*std_im_pref_2W10
  wim_prefW11 =~ 1*std_im_pref_2W11
  wim_prefW14 =~ 1*std_im_pref_2W14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 + wim_prefW10 ~ wstd_auth01W7 + wim_prefW7
  wstd_auth01W11 + wim_prefW11 ~ wstd_auth01W10 + wim_prefW10
  wstd_auth01W14 + wim_prefW14 ~ wstd_auth01W11 + wim_prefW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wim_prefW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wim_prefW10
  wstd_auth01W11 ~~ wim_prefW11
  wstd_auth01W14 ~~ wim_prefW14
  
  
  # Constrain the variance and covariance of the random intercepts to 0 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIim_prefW ~~ 0*RIim_prefW
  RIstd_auth01W ~~ 0*RIim_prefW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wim_prefW7 ~~ wim_prefW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wim_prefW10 ~~ wim_prefW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wim_prefW11 ~~ wim_prefW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wim_prefW14 ~~ wim_prefW14 
  
# Fix error variance of indicators to (1 - alpha)*var
std_auth01W7 ~~ 0.02*std_auth01W7
std_auth01W10 ~~ 0.016*std_auth01W10
std_auth01W11 ~~ 0.018*std_auth01W11
std_auth01W14 ~~ 0.019*std_auth01W14
std_im_pref_2W7 ~~ 0.0147*std_im_pref_2W7
std_im_pref_2W10 ~~ 0.000155*std_im_pref_2W10
std_im_pref_2W11 ~~ 0.000147*std_im_pref_2W11
std_im_pref_2W14 ~~ 0.000512*std_im_pref_2W14

  
  # incorporate tvc - age
wstd_auth01W7 + wim_prefW7 ~ ageW7
wstd_auth01W10 + wim_prefW10 ~ ageW10
wstd_auth01W11 + wim_prefW11 ~ ageW11
wstd_auth01W14 + wim_prefW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc - ethnicity
wstd_auth01W7 + wim_prefW7 ~ ethnic7
wstd_auth01W10 + wim_prefW10 ~ ethnic10
wstd_auth01W11 + wim_prefW11 ~ ethnic11
wstd_auth01W14 + wim_prefW14 ~ ethnic14

ethnic7 ~~ ethnic7
ethnic10 ~~ ethnic10
ethnic11 ~~ ethnic11
ethnic14 ~~ ethnic14

ethnic7 ~~ ethnic10 + ethnic11 + ethnic14
ethnic10 ~~ ethnic11 + ethnic14
ethnic11 ~~ ethnic14

# incorporate tvc – social grade
wstd_auth01W7 + wim_prefW7 ~ p_socgradeW7
wstd_auth01W10 + wim_prefW10 ~ p_socgradeW10
wstd_auth01W11 + wim_prefW11 ~ p_socgradeW11
wstd_auth01W14 + wim_prefW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – education
wstd_auth01W7 + wim_prefW7 ~ higheduc7
wstd_auth01W10 + wim_prefW10 ~ higheduc10
wstd_auth01W11 + wim_prefW11 ~ higheduc11
wstd_auth01W14 + wim_prefW14 ~ higheduc14

higheduc7 ~~ higheduc7
higheduc10 ~~ higheduc10
higheduc11 ~~ higheduc11
higheduc14 ~~ higheduc14

higheduc7 ~~ higheduc10 + higheduc11 + higheduc14
higheduc10 ~~ higheduc11 + higheduc14
higheduc11 ~~ higheduc14

# incorporate tvc – income
wstd_auth01W7 + wim_prefW7 ~ incomerec7
wstd_auth01W10 + wim_prefW10 ~ incomerec10
wstd_auth01W11 + wim_prefW11 ~ incomerec11
wstd_auth01W14 + wim_prefW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14
  '
std_immig_CLPM_L_noEQ.fit <- lavaan(std_immig_CLPM_L_noEQ, 
                                        data = att_dat_low,
                                        estimator = "MLR",
                                        missing = 'ML',
                                        meanstructure = T, 
                                        int.ov.free = T) 
summary(std_immig_CLPM_L_noEQ.fit, standardized = T)

# Model fit 

fitmeasures(std_immig_CLPM_H.fit)
fitmeasures(std_immig_CLPM_H_noEQ.fit)

# Model fit for high engaged, constrained. 
# AIC = 183907.886; CFI = 0.971; RMSEA = 0.043; SRMR = 0.087.
# Model fit for high engaged, unconstrained. 
# AIC = 183862.004; CFI = 0.971; RMSEA = 0.043; SRMR = 0.087.

fitmeasures(std_immig_CLPM_L.fit)
fitmeasures(std_immig_CLPM_L_noEQ.fit)

# Model fit for less engaged, constrained. 
# AIC = 138691.011; CFI = 0.98; RMSEA = 0.033; SRMR = 0.095.
# Model fit for less engaged, unconstrained. 
# AIC = 138671.369; CFI = 0.98; RMSEA = 0.033; SRMR = 0.094.

# Results

standardizedSolution(std_immig_CLPM_H.fit)
standardizedSolution(std_immig_CLPM_H_noEQ.fit)

standardizedSolution(std_immig_CLPM_L.fit)
standardizedSolution(std_immig_CLPM_L_noEQ.fit)



#### 16. Multiple group CLPMs - EU -----------

# First - fit a model with unconstrained differences between groups. 

std_eu_attU_CLPM <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wstd_eu_intW7 =~ 1*std_eu_intW7
  wstd_eu_intW10 =~ 1*std_eu_intW10
  wstd_eu_intW11 =~ 1*std_eu_intW11
  wstd_eu_intW14 =~ 1*std_eu_intW14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 + wstd_eu_intW10 ~ wstd_auth01W7 + wstd_eu_intW7
  wstd_auth01W11 + wstd_eu_intW11 ~ wstd_auth01W10 + wstd_eu_intW10
  wstd_auth01W14 + wstd_eu_intW14 ~ wstd_auth01W11 + wstd_eu_intW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wstd_eu_intW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wstd_eu_intW10
  wstd_auth01W11 ~~ wstd_eu_intW11
  wstd_auth01W14 ~~ wstd_eu_intW14
  
  
  # Fix the variance and covariance of the random intercepts to 0 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIstd_eu_intW ~~ 0*RIstd_eu_intW
  RIstd_auth01W ~~ 0*RIstd_eu_intW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ wstd_eu_intW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ wstd_eu_intW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wstd_eu_intW11 ~~ wstd_eu_intW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wstd_eu_intW14 ~~ wstd_eu_intW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W10 ~~ 0.0058*std_auth01W10
std_auth01W11 ~~ 0.002*std_auth01W11
std_auth01W14 ~~ 0.00725*std_auth01W14
std_eu_intW10 ~~ 0.0028*std_eu_intW10
std_eu_intW11 ~~ 0.00194*std_eu_intW11
std_eu_intW14 ~~ 0*std_eu_intW14


'

std_eu_attU_CLPM.fit <- lavaan(std_eu_attU_CLPM, 
                               estimator = "MLR",
                               data = att_dat,
                               group = "polatt_group",
                               missing = 'ML', 
                               meanstructure = T, 
                               int.ov.free = T) 
summary(std_eu_attU_CLPM.fit, standardized = T)

# Then a model with constrained differences between groups

std_eu_attC_CLPM <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wstd_eu_intW7 =~ 1*std_eu_intW7
  wstd_eu_intW10 =~ 1*std_eu_intW10
  wstd_eu_intW11 =~ 1*std_eu_intW11
  wstd_eu_intW14 =~ 1*std_eu_intW14
  
  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 ~ c(a1, a1)*wstd_auth01W7 + c(b1, b1)*wstd_eu_intW7
  wstd_eu_intW10 ~ c(c1, c1)*wstd_auth01W7 + c(d1, d1)*wstd_eu_intW7
  wstd_auth01W11 ~ c(a2, a2)*wstd_auth01W10 + c(b2, b2)*wstd_eu_intW10
  wstd_eu_intW11 ~ c(c2, c2)*wstd_auth01W10 + c(d2, d2)*wstd_eu_intW10
  wstd_auth01W14 ~ c(a3, a3)*wstd_auth01W11 + c(b3, b3)*wstd_eu_intW11
  wstd_eu_intW14 ~ c(c3, c3)*wstd_auth01W11 + c(d3, d3)*wstd_eu_intW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wstd_eu_intW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wstd_eu_intW10
  wstd_auth01W11 ~~ wstd_eu_intW11
  wstd_auth01W14 ~~ wstd_eu_intW14
  
  # Fix the variance and covariance of the random intercepts to 0 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIstd_eu_intW ~~ 0*RIstd_eu_intW
  RIstd_auth01W ~~ 0*RIstd_eu_intW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ wstd_eu_intW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ wstd_eu_intW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wstd_eu_intW11 ~~ wstd_eu_intW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wstd_eu_intW14 ~~ wstd_eu_intW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W10 ~~ 0.0058*std_auth01W10
std_auth01W11 ~~ 0.002*std_auth01W11
std_auth01W14 ~~ 0.00725*std_auth01W14
std_eu_intW10 ~~ 0.0028*std_eu_intW10
std_eu_intW11 ~~ 0.00194*std_eu_intW11
std_eu_intW14 ~~ 0*std_eu_intW14
'

std_eu_attC_CLPM.fit <- lavaan(std_eu_attC_CLPM, 
                               estimator = "MLR",
                               data = att_dat,
                               group = "polatt_group",
                               missing = 'ML', 
                               meanstructure = T, 
                               int.ov.free = T) 
summary(std_eu_attC_CLPM.fit, standardized = T)

# ANOVA to compare model fit 

anova(std_eu_attC_CLPM.fit, std_eu_attU_CLPM.fit)

# p < 0.05, so the high and low engagement groups do differ. Now I fit separate 
# models for the two groups so that equality constraints on the parameters
# over time can be imposed 

## Highly engaged group 

# with equality constraints 

std_eu_CLPM_H <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ NA*std_auth01W7
  wstd_auth01W10 =~ NA*std_auth01W10
  wstd_auth01W11 =~ NA*std_auth01W11 
  wstd_auth01W14 =~ NA*std_auth01W14

  wstd_eu_intW7 =~ NA*std_eu_intW7
  wstd_eu_intW10 =~ NA*std_eu_intW10
  wstd_eu_intW11 =~ NA*std_eu_intW11
  wstd_eu_intW14 =~ NA*std_eu_intW14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 ~ a*wstd_auth01W7 + b*wstd_eu_intW7
  wstd_eu_intW10 ~ c*wstd_auth01W7 + d*wstd_eu_intW7
  wstd_auth01W11 ~ a*wstd_auth01W10 + b*wstd_eu_intW10
  wstd_eu_intW11 ~ c*wstd_auth01W10 + d*wstd_eu_intW10
  wstd_auth01W14 ~ a*wstd_auth01W11 + b*wstd_eu_intW11
  wstd_eu_intW14 ~ c*wstd_auth01W11 + d*wstd_eu_intW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ cor1*wstd_eu_intW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ rcov2*wstd_eu_intW10
  wstd_auth01W11 ~~ rcov3*wstd_eu_intW11
  wstd_auth01W14 ~~ rcov4*wstd_eu_intW14
  
  
  # Fix the variance and covariance of the random intercepts to 0 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIstd_eu_intW ~~ 0*RIstd_eu_intW
  RIstd_auth01W ~~ 0*RIstd_eu_intW

  # Set variances of within-components at first wave to 1
  wstd_auth01W7 ~~ 1*wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ 1*wstd_eu_intW7 
  
  wstd_auth01W10 ~~ rvx2*wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ rvy2*wstd_eu_intW10 
  wstd_auth01W11 ~~ rvx3*wstd_auth01W11
  wstd_eu_intW11 ~~ rvy3*wstd_eu_intW11 
  wstd_auth01W14 ~~ rvx4*wstd_auth01W14 
  wstd_eu_intW14 ~~ rvy4*wstd_eu_intW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W10 ~~ 0.0058*std_auth01W10
std_auth01W11 ~~ 0.002*std_auth01W11
std_auth01W14 ~~ 0.00725*std_auth01W14
std_eu_intW10 ~~ 0.0028*std_eu_intW10
std_eu_intW11 ~~ 0.00194*std_eu_intW11
std_eu_intW14 ~~ 0*std_eu_intW14

 # Compute correlations of within-components at each wave
  cor2 := a*c + b*d + a*d*cor1 + b*c*cor1 + rcov2
  cor3 := a*c + b*d + a*d*cor2 + b*c*cor2 + rcov3
  cor4 := a*c + b*d + a*d*cor3 + b*c*cor3 + rcov4
 
  # Contrain residual variances of within-components such that variance of each 
  # within-component equals 1
  rvx2 == 1 - (a*a + b*b + 2*a*b*cor1)
  rvy2 == 1 - (c*c + d*d + 2*c*d*cor1)
  rvx3 == 1 - (a*a + b*b + 2*a*b*cor2)
  rvy3 == 1 - (c*c + d*d + 2*c*d*cor2)
  rvx4 == 1 - (a*a + b*b + 2*a*b*cor3)
  rvy4 == 1 - (c*c + d*d + 2*c*d*cor3)
  
  # incorporate tvc - age
wstd_auth01W7 + wstd_eu_intW7 ~ ageW7
wstd_auth01W10 + wstd_eu_intW10 ~ ageW10
wstd_auth01W11 + wstd_eu_intW11 ~ ageW11
wstd_auth01W14 + wstd_eu_intW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc - ethnicity
wstd_auth01W7 + wstd_eu_intW7 ~ ethnic7
wstd_auth01W10 + wstd_eu_intW10 ~ ethnic10
wstd_auth01W11 + wstd_eu_intW11 ~ ethnic11
wstd_auth01W14 + wstd_eu_intW14 ~ ethnic14

ethnic7 ~~ ethnic7
ethnic10 ~~ ethnic10
ethnic11 ~~ ethnic11
ethnic14 ~~ ethnic14

ethnic7 ~~ ethnic10 + ethnic11 + ethnic14
ethnic10 ~~ ethnic11 + ethnic14
ethnic11 ~~ ethnic14

# incorporate tvc – social grade
wstd_auth01W7 + wstd_eu_intW7 ~ p_socgradeW7
wstd_auth01W10 + wstd_eu_intW10 ~ p_socgradeW10
wstd_auth01W11 + wstd_eu_intW11 ~ p_socgradeW11
wstd_auth01W14 + wstd_eu_intW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – education
wstd_auth01W7 + wstd_eu_intW7 ~ higheduc7
wstd_auth01W10 + wstd_eu_intW10 ~ higheduc10
wstd_auth01W11 + wstd_eu_intW11 ~ higheduc11
wstd_auth01W14 + wstd_eu_intW14 ~ higheduc14

higheduc7 ~~ higheduc7
higheduc10 ~~ higheduc10
higheduc11 ~~ higheduc11
higheduc14 ~~ higheduc14

higheduc7 ~~ higheduc10 + higheduc11 + higheduc14
higheduc10 ~~ higheduc11 + higheduc14
higheduc11 ~~ higheduc14

# incorporate tvc – income
wstd_auth01W7 + wstd_eu_intW7 ~ incomerec7
wstd_auth01W10 + wstd_eu_intW10 ~ incomerec10
wstd_auth01W11 + wstd_eu_intW11 ~ incomerec11
wstd_auth01W14 + wstd_eu_intW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14
'

std_eu_CLPM_H.fit <- lavaan(std_eu_CLPM_H, 
                                estimator = "MLR",
                                data = att_dat_high,
                                missing = 'ML', 
                                meanstructure = T, 
                                int.ov.free = T) 
summary(std_eu_CLPM_H.fit, standardized = T)


# Without equality constraints 

std_eu_CLPM_H_noEQ <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wstd_eu_intW7 =~ 1*std_eu_intW7
  wstd_eu_intW10 =~ 1*std_eu_intW10
  wstd_eu_intW11 =~ 1*std_eu_intW11
  wstd_eu_intW14 =~ 1*std_eu_intW14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 + wstd_eu_intW10 ~ wstd_auth01W7 + wstd_eu_intW7
  wstd_auth01W11 + wstd_eu_intW11 ~ wstd_auth01W10 + wstd_eu_intW10
  wstd_auth01W14 + wstd_eu_intW14 ~ wstd_auth01W11 + wstd_eu_intW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wstd_eu_intW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wstd_eu_intW10
  wstd_auth01W11 ~~ wstd_eu_intW11
  wstd_auth01W14 ~~ wstd_eu_intW14
  
  
  # Fix the variance and covariance of the random intercepts to 0 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIstd_eu_intW ~~ 0*RIstd_eu_intW
  RIstd_auth01W ~~ 0*RIstd_eu_intW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ wstd_eu_intW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ wstd_eu_intW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wstd_eu_intW11 ~~ wstd_eu_intW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wstd_eu_intW14 ~~ wstd_eu_intW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W10 ~~ 0.0058*std_auth01W10
std_auth01W11 ~~ 0.002*std_auth01W11
std_auth01W14 ~~ 0.00725*std_auth01W14
std_eu_intW10 ~~ 0.0028*std_eu_intW10
std_eu_intW11 ~~ 0.00194*std_eu_intW11
std_eu_intW14 ~~ 0*std_eu_intW14

  # incorporate tvc - age
wstd_auth01W7 + wstd_eu_intW7 ~ ageW7
wstd_auth01W10 + wstd_eu_intW10 ~ ageW10
wstd_auth01W11 + wstd_eu_intW11 ~ ageW11
wstd_auth01W14 + wstd_eu_intW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc - ethnicity
wstd_auth01W7 + wstd_eu_intW7 ~ ethnic7
wstd_auth01W10 + wstd_eu_intW10 ~ ethnic10
wstd_auth01W11 + wstd_eu_intW11 ~ ethnic11
wstd_auth01W14 + wstd_eu_intW14 ~ ethnic14

ethnic7 ~~ ethnic7
ethnic10 ~~ ethnic10
ethnic11 ~~ ethnic11
ethnic14 ~~ ethnic14

ethnic7 ~~ ethnic10 + ethnic11 + ethnic14
ethnic10 ~~ ethnic11 + ethnic14
ethnic11 ~~ ethnic14

# incorporate tvc – social grade
wstd_auth01W7 + wstd_eu_intW7 ~ p_socgradeW7
wstd_auth01W10 + wstd_eu_intW10 ~ p_socgradeW10
wstd_auth01W11 + wstd_eu_intW11 ~ p_socgradeW11
wstd_auth01W14 + wstd_eu_intW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – education
wstd_auth01W7 + wstd_eu_intW7 ~ higheduc7
wstd_auth01W10 + wstd_eu_intW10 ~ higheduc10
wstd_auth01W11 + wstd_eu_intW11 ~ higheduc11
wstd_auth01W14 + wstd_eu_intW14 ~ higheduc14

higheduc7 ~~ higheduc7
higheduc10 ~~ higheduc10
higheduc11 ~~ higheduc11
higheduc14 ~~ higheduc14

higheduc7 ~~ higheduc10 + higheduc11 + higheduc14
higheduc10 ~~ higheduc11 + higheduc14
higheduc11 ~~ higheduc14

# incorporate tvc – income
wstd_auth01W7 + wstd_eu_intW7 ~ incomerec7
wstd_auth01W10 + wstd_eu_intW10 ~ incomerec10
wstd_auth01W11 + wstd_eu_intW11 ~ incomerec11
wstd_auth01W14 + wstd_eu_intW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14
'

std_eu_CLPM_H_noEQ.fit <- lavaan(std_eu_CLPM_H_noEQ, 
                                     estimator = "MLR",
                                     data = att_dat_high,
                                     missing = 'ML', 
                                     meanstructure = T, 
                                     int.ov.free = T) 
summary(std_eu_CLPM_H_noEQ.fit, standardized = T)


## Less engaged group 

# With equality constraints 

std_eu_CLPM_L <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ NA*std_auth01W7
  wstd_auth01W10 =~ NA*std_auth01W10
  wstd_auth01W11 =~ NA*std_auth01W11 
  wstd_auth01W14 =~ NA*std_auth01W14

  wstd_eu_intW7 =~ NA*std_eu_intW7
  wstd_eu_intW10 =~ NA*std_eu_intW10
  wstd_eu_intW11 =~ NA*std_eu_intW11
  wstd_eu_intW14 =~ NA*std_eu_intW14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 ~ a*wstd_auth01W7 + b*wstd_eu_intW7
  wstd_eu_intW10 ~ c*wstd_auth01W7 + d*wstd_eu_intW7
  wstd_auth01W11 ~ a*wstd_auth01W10 + b*wstd_eu_intW10
  wstd_eu_intW11 ~ c*wstd_auth01W10 + d*wstd_eu_intW10
  wstd_auth01W14 ~ a*wstd_auth01W11 + b*wstd_eu_intW11
  wstd_eu_intW14 ~ c*wstd_auth01W11 + d*wstd_eu_intW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ cor1*wstd_eu_intW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ rcov2*wstd_eu_intW10
  wstd_auth01W11 ~~ rcov3*wstd_eu_intW11
  wstd_auth01W14 ~~ rcov4*wstd_eu_intW14
  
  
  # Fix the variance and covariance of the random intercepts to 0 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIstd_eu_intW ~~ 0*RIstd_eu_intW
  RIstd_auth01W ~~ 0*RIstd_eu_intW

  # Set variances of within-components at first wave to 1
  wstd_auth01W7 ~~ 1*wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ 1*wstd_eu_intW7 
  
  wstd_auth01W10 ~~ rvx2*wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ rvy2*wstd_eu_intW10 
  wstd_auth01W11 ~~ rvx3*wstd_auth01W11
  wstd_eu_intW11 ~~ rvy3*wstd_eu_intW11 
  wstd_auth01W14 ~~ rvx4*wstd_auth01W14 
  wstd_eu_intW14 ~~ rvy4*wstd_eu_intW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W10 ~~ 0.0058*std_auth01W10
std_auth01W11 ~~ 0.002*std_auth01W11
std_auth01W14 ~~ 0.00725*std_auth01W14
std_eu_intW10 ~~ 0.0028*std_eu_intW10
std_eu_intW11 ~~ 0.00194*std_eu_intW11
std_eu_intW14 ~~ 0*std_eu_intW14

 # Compute correlations of within-components at each wave
  cor2 := a*c + b*d + a*d*cor1 + b*c*cor1 + rcov2
  cor3 := a*c + b*d + a*d*cor2 + b*c*cor2 + rcov3
  cor4 := a*c + b*d + a*d*cor3 + b*c*cor3 + rcov4
 
  # Contrain residual variances of within-components such that variance of each 
  # within-component equals 1
  rvx2 == 1 - (a*a + b*b + 2*a*b*cor1)
  rvy2 == 1 - (c*c + d*d + 2*c*d*cor1)
  rvx3 == 1 - (a*a + b*b + 2*a*b*cor2)
  rvy3 == 1 - (c*c + d*d + 2*c*d*cor2)
  rvx4 == 1 - (a*a + b*b + 2*a*b*cor3)
  rvy4 == 1 - (c*c + d*d + 2*c*d*cor3)
  
    # incorporate tvc - age
wstd_auth01W7 + wstd_eu_intW7 ~ ageW7
wstd_auth01W10 + wstd_eu_intW10 ~ ageW10
wstd_auth01W11 + wstd_eu_intW11 ~ ageW11
wstd_auth01W14 + wstd_eu_intW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc - ethnicity
wstd_auth01W7 + wstd_eu_intW7 ~ ethnic7
wstd_auth01W10 + wstd_eu_intW10 ~ ethnic10
wstd_auth01W11 + wstd_eu_intW11 ~ ethnic11
wstd_auth01W14 + wstd_eu_intW14 ~ ethnic14

ethnic7 ~~ ethnic7
ethnic10 ~~ ethnic10
ethnic11 ~~ ethnic11
ethnic14 ~~ ethnic14

ethnic7 ~~ ethnic10 + ethnic11 + ethnic14
ethnic10 ~~ ethnic11 + ethnic14
ethnic11 ~~ ethnic14

# incorporate tvc – social grade
wstd_auth01W7 + wstd_eu_intW7 ~ p_socgradeW7
wstd_auth01W10 + wstd_eu_intW10 ~ p_socgradeW10
wstd_auth01W11 + wstd_eu_intW11 ~ p_socgradeW11
wstd_auth01W14 + wstd_eu_intW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – education
wstd_auth01W7 + wstd_eu_intW7 ~ higheduc7
wstd_auth01W10 + wstd_eu_intW10 ~ higheduc10
wstd_auth01W11 + wstd_eu_intW11 ~ higheduc11
wstd_auth01W14 + wstd_eu_intW14 ~ higheduc14

higheduc7 ~~ higheduc7
higheduc10 ~~ higheduc10
higheduc11 ~~ higheduc11
higheduc14 ~~ higheduc14

higheduc7 ~~ higheduc10 + higheduc11 + higheduc14
higheduc10 ~~ higheduc11 + higheduc14
higheduc11 ~~ higheduc14

# incorporate tvc – income
wstd_auth01W7 + wstd_eu_intW7 ~ incomerec7
wstd_auth01W10 + wstd_eu_intW10 ~ incomerec10
wstd_auth01W11 + wstd_eu_intW11 ~ incomerec11
wstd_auth01W14 + wstd_eu_intW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14
'

std_eu_CLPM_L.fit <- lavaan(std_eu_CLPM_L, 
                                estimator = "MLR",
                                data = att_dat_low,
                                missing = 'ML', 
                                meanstructure = T, 
                                int.ov.free = T) 
summary(std_eu_CLPM_L.fit, standardized = T)

# Without equality constraints 

std_eu_CLPM_L_noEQ <- '
  # Create between components (random intercepts)
  RIstd_auth01W  =~ 1*std_auth01W7 + 1*std_auth01W10 + 1*std_auth01W11 + 1*std_auth01W14 
  RIstd_eu_intW =~ 1*std_eu_intW7 + 1*std_eu_intW10 + 1*std_eu_intW11 + 1*std_eu_intW14 
  
  # Create within-person centered variables
  wstd_auth01W7 =~ 1*std_auth01W7
  wstd_auth01W10 =~ 1*std_auth01W10
  wstd_auth01W11 =~ 1*std_auth01W11 
  wstd_auth01W14 =~ 1*std_auth01W14

  wstd_eu_intW7 =~ 1*std_eu_intW7
  wstd_eu_intW10 =~ 1*std_eu_intW10
  wstd_eu_intW11 =~ 1*std_eu_intW11
  wstd_eu_intW14 =~ 1*std_eu_intW14

  # Estimate the lagged effects between the within-person centered variables.
  wstd_auth01W10 + wstd_eu_intW10 ~ wstd_auth01W7 + wstd_eu_intW7
  wstd_auth01W11 + wstd_eu_intW11 ~ wstd_auth01W10 + wstd_eu_intW10
  wstd_auth01W14 + wstd_eu_intW14 ~ wstd_auth01W11 + wstd_eu_intW11

  # Estimate the covariance between the within-person centered variables at the first wave. 
  wstd_auth01W7 ~~ wstd_eu_intW7 # Covariance
  
  # Estimate the covariances between the residuals of the within-person centered variables (the innovations).
  wstd_auth01W10 ~~ wstd_eu_intW10
  wstd_auth01W11 ~~ wstd_eu_intW11
  wstd_auth01W14 ~~ wstd_eu_intW14
  
  
  # Fix the variance and covariance of the random intercepts to 0 
  RIstd_auth01W ~~ 0*RIstd_auth01W
  RIstd_eu_intW ~~ 0*RIstd_eu_intW
  RIstd_auth01W ~~ 0*RIstd_eu_intW

  # Estimate the (residual) variance of the within-person centered variables.
  wstd_auth01W7 ~~ wstd_auth01W7 # Variances
  wstd_eu_intW7 ~~ wstd_eu_intW7 
  wstd_auth01W10 ~~ wstd_auth01W10 # Residual variances
  wstd_eu_intW10 ~~ wstd_eu_intW10 
  wstd_auth01W11 ~~ wstd_auth01W11
  wstd_eu_intW11 ~~ wstd_eu_intW11 
  wstd_auth01W14 ~~ wstd_auth01W14 
  wstd_eu_intW14 ~~ wstd_eu_intW14 
  
  # Fix error variance of indicators to (1 - alpha)*var
std_auth01W10 ~~ 0.0058*std_auth01W10
std_auth01W11 ~~ 0.002*std_auth01W11
std_auth01W14 ~~ 0.00725*std_auth01W14
std_eu_intW10 ~~ 0.0028*std_eu_intW10
std_eu_intW11 ~~ 0.00194*std_eu_intW11
std_eu_intW14 ~~ 0*std_eu_intW14

  # incorporate tvc - age
wstd_auth01W7 + wstd_eu_intW7 ~ ageW7
wstd_auth01W10 + wstd_eu_intW10 ~ ageW10
wstd_auth01W11 + wstd_eu_intW11 ~ ageW11
wstd_auth01W14 + wstd_eu_intW14 ~ ageW14

ageW7 ~~ ageW7
ageW10 ~~ ageW10
ageW11 ~~ ageW11
ageW14 ~~ ageW14

ageW7 ~~ ageW10 + ageW11 + ageW14
ageW10 ~~ ageW11 + ageW14
ageW11 ~~ ageW14

# incorporate tvc - ethnicity
wstd_auth01W7 + wstd_eu_intW7 ~ ethnic7
wstd_auth01W10 + wstd_eu_intW10 ~ ethnic10
wstd_auth01W11 + wstd_eu_intW11 ~ ethnic11
wstd_auth01W14 + wstd_eu_intW14 ~ ethnic14

ethnic7 ~~ ethnic7
ethnic10 ~~ ethnic10
ethnic11 ~~ ethnic11
ethnic14 ~~ ethnic14

ethnic7 ~~ ethnic10 + ethnic11 + ethnic14
ethnic10 ~~ ethnic11 + ethnic14
ethnic11 ~~ ethnic14

# incorporate tvc – social grade
wstd_auth01W7 + wstd_eu_intW7 ~ p_socgradeW7
wstd_auth01W10 + wstd_eu_intW10 ~ p_socgradeW10
wstd_auth01W11 + wstd_eu_intW11 ~ p_socgradeW11
wstd_auth01W14 + wstd_eu_intW14 ~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW7
p_socgradeW10 ~~ p_socgradeW10
p_socgradeW11 ~~ p_socgradeW11
p_socgradeW14 ~~ p_socgradeW14

p_socgradeW7 ~~ p_socgradeW10 + p_socgradeW11 + p_socgradeW14
p_socgradeW10 ~~ p_socgradeW11 + p_socgradeW14
p_socgradeW11 ~~ p_socgradeW14

# incorporate tvc – education
wstd_auth01W7 + wstd_eu_intW7 ~ higheduc7
wstd_auth01W10 + wstd_eu_intW10 ~ higheduc10
wstd_auth01W11 + wstd_eu_intW11 ~ higheduc11
wstd_auth01W14 + wstd_eu_intW14 ~ higheduc14

higheduc7 ~~ higheduc7
higheduc10 ~~ higheduc10
higheduc11 ~~ higheduc11
higheduc14 ~~ higheduc14

higheduc7 ~~ higheduc10 + higheduc11 + higheduc14
higheduc10 ~~ higheduc11 + higheduc14
higheduc11 ~~ higheduc14

# incorporate tvc – income
wstd_auth01W7 + wstd_eu_intW7 ~ incomerec7
wstd_auth01W10 + wstd_eu_intW10 ~ incomerec10
wstd_auth01W11 + wstd_eu_intW11 ~ incomerec11
wstd_auth01W14 + wstd_eu_intW14 ~ incomerec14

incomerec7 ~~ incomerec7
incomerec10 ~~ incomerec10
incomerec11 ~~ incomerec11
incomerec14 ~~ incomerec14

incomerec7 ~~ incomerec10 + incomerec11 + incomerec14
incomerec10 ~~ incomerec11 + incomerec14
incomerec11 ~~ incomerec14

'

std_eu_CLPM_L_noEQ.fit <- lavaan(std_eu_CLPM_L_noEQ, 
                                     estimator = "MLR",
                                     data = att_dat_low,
                                     missing = 'ML', 
                                     meanstructure = T, 
                                     int.ov.free = T) 
summary(std_eu_CLPM_L_noEQ.fit, standardized = T)

# Model fit 

fitmeasures(std_eu_CLPM_H.fit)
fitmeasures(std_eu_CLPM_H_noEQ.fit)

# Model fit for high engaged, constrained. 
# AIC = 185568.195; CFI = 0.97; RMSEA = 0.043; SRMR = 0.088.
# Model fit for high engaged, unconstrained. 
# AIC = 185432.045; CFI = 0.971; RMSEA = 0.043; SRMR = 0.087.

fitmeasures(std_eu_CLPM_L.fit)
fitmeasures(std_eu_CLPM_L_noEQ.fit)

# Model fit for less engaged, constrained. 
# AIC = 141302.346; CFI = 0.979; RMSEA = 0.034; SRMR = 0.094.  
# Model fit for less engaged, unconstrained. 
# AIC = 141245.976; CFI = 0.98; RMSEA = 0.034; SRMR = 0.094.

# Results

standardizedSolution(std_eu_CLPM_H.fit)
standardizedSolution(std_eu_CLPM_H_noEQ.fit)

standardizedSolution(std_eu_CLPM_L.fit)
standardizedSolution(std_eu_CLPM_L_noEQ.fit)